DIAGNOSTIC SIGNAL PROCESSING METHOD AND SYSTEM 

RELATED APPLICATION INFORMATION 

This application claims the benefit of U.S. Provisional Patent Application Serial No. 
60/426,348, filed on November 14, 2002; 

FIELD OF THE INVENTION 

The present invention relates generally to signal processing applications and, in 
particular, to signal processing applications involving time-to-frequency domain transforms 
5 where it is useful to define a spectrum including frequency information for frequencies other than 
a fimdamental frequency and harmonies thereof, or where it is useful to generate a spectrum 
based on a portion of the time-based signal representing less than a period of a signal component. 
The invention has particular application to methods wherein a received signal is decomposed or 
otherwise processed in a manner that provides information regarding properties of a material 
10 modulating or otherwise producing the received signal. The received signal may emanate firom 
tissue interacting with an interrogating signal that was introduced into an organism and was 
modified from the interrogating signal by one or more influences of the organism's material, e.g., 
organs, tissue, fluids or other structure. 

15 BACKGROUND OF THE INVENTION 

Noninvasive diagnoses using interrogating signals are frequently employed to ascertain 
the condition of tissue and functional parameters of organs. Such interrogating signals may use, 
for example, electromagnetic or sonic, including ultrasonic, energy. One example is transmitting 
ultrasonic waves into an organism, such as a human being, and receiving a retum signal that is a 
20 modified version of the transmitted interrogating signal, with the modifications caused by the 
properties of one or more tissue types. The received signal is measured and processed to 
evaluate the condition or properties of the tissue(s) or other material(s) with which the 
interrogating signal interacts. In particular, interrogating ultrasonic signals are often transmitted 
into an organism and the received signal is due to backscattering or echos of the interrogating 
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signal by one or more tissue types. Even more particularly, an ultrasonic interrogating signal is 
often transmitted into a region of the organism where blood flow is expected and the received 
signal is a modified version of the interrogating signal wherein such modifications are due in 
whole or part to frequency shifts caused by Doppler effects caused by movement of tissues, 
5 including movement of blood. In this case it is desirable to measure and process the received 

signal in a manner such that the amount of Doppler frequency shift that occurs can be ascertained 
in order to estimate the various velocities of the moving tissue, particularly the velocities of the 
flowing blood. 

Ultrasonic Doppler techniques typically entail using a hand-held probe that contains one 

10 or more transducers that are used to transmit the interrogating signal and to produce a response, 
typically in the forrri of a time-varying voltage, to the received signal. The interrogating 
ultrasonic signal typically has a frequency in the radio frequency (RF) range of about 0.5 million 
^ Hertz (0.5 MHz) to about 20 MHz, although other frequencies may be used. The interrogating 
signal may be sent continuously, in which case the method employed is termed "continuous 

15 wave" (CW) or the interrogating signal may be sent intermittently, usually on a periodic basis, in 
which case the method employed is termed "pulsed." When pulsed Doppler is used it is possible 
to interrogate tissue over a range of specified distances from the transducers by measuring the 
received signal only for a selected duration that occurs after a predetermined time delay after the 
pulse is sent. If the speed of sound for the intervening tissue is known then the time delay 

20 converts directly into a depth of tissue being interrogated. Interrogating a specified tissue depth 
often provides clinically valuable information that cannot be obtained if CW Doppler is used. 
CW Doppler systems are often simpler than pulsed Doppler systems, although they cannot 
provide information that is known to come from tissue associated with a specific time delay (and 
corresponding estimated tissue depth). 

25 Regardless of whether the Doppler system is pulsed or continuous, the spectral content of 

the received signal is analyzed to assess how the tissue modified the interrogating signal. The 
spectral frequency content of the signal is a measure of the amplitudes of each of the constituent 
frequencies that, when combined, form the total signal. The spectral content is typically 
represented as either the amplitude of each of the frequencies or as the po>yer of each of the 
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frequencies. When the spectral content data are graphed the resulting graphs have a variety of 
names including spectrograms and spectral density functions. The term "spectral density 
function" (SD function) will be used as the general term below. When the frequency content of 
the interrogating signal is known it can be compared to the spectrum of the received signal and 
the differences between the spectra can be interpreted as arising from the velocity profiles of the 
tissues through which the signals passed. The velocity profiles can be back-calculated from 
spectral differences using established methods that are well known. 

To date, a variety of methods have been used to estimate the spectral frequency content 
of the received signals. The main spectral analysis methods are traditionally classified as being 
either nonparametric or parametric methods. 

The nonparametric methods aire related to using Fourier series and they estimate the 
spectral content at equally spaced fi^equency points. The extant methods include Fourier and 
z-transform (including the Chirp z-transform). For brevity, all methods which employ assessing 
spectral content at equally spaced frequency points are referred to as Fourier methods. 

Fourier methods, typically in the form of fast Fourier transforms (FFT), are well known 
and decompose a signal into frequencies that are integral multiples of a fundamental frequency. 
The integral multiple fi-equencies form an orthogonal basis in the frequency domain. This 
information can be used to approximate the signal as a trigonometric polynomial, as is well 
known. Fourier methods have been widely used to analyze the received signal because such 
methods have been refined to the point where they are fast and have well defined properties. 

The other class of spectral analysis methods are parametric methods. Parametric methods 
use time-series analysis to estimate the paranieters for a rational function model. The rational 
function can have either only poles,. which is an autoregressive (AR) model, only zeros, which is 
a moving average (MA) model, or both poles and zeros, which is an autoregressive-moving 
average (ARMA) model. Included in these methods are maximum entropy methods (MEM). The 
number of model parameters used in a parametric method is termed the order of the model. The 
various methods that exist for calculating the values of the model parameters, such as the Burg 
algorithm and Durbin's first and second methods, are well known. An advantage of parametric 
methods is that they can produce reasonable spectral estimates using very few cycles or fi'actions 



of cycles of data. 

SUMMARY OF THE INVENTION 
It has been recognized that alternate spectral analysis methods including alternate time-to^ 
5 . frequency domain transforms may be beneficial for a variety of signal processing applications, 
including certain medical applications. In particular, in spite of the extensive use of Fourier 
methods, they suffer from significant deficiencies. Parametric methods also suffer significant 
deficiencies. With regard to Fourier methods, the most significant deficiency is that these 
methods produce spectral estimates for only a small number of discrete frequencies, even when it 

10 is expected or known that the real spectrum covers a continuous range of frequencies or that the 
spectrum may contain frequencies different from those that form the basis that is the output of 
the Fourier transform. For example, the signal may contain frequencies that are not an integral 
multiple of the fiindamental frequency. When the basis is not correct, either because the 
underlying spectrum is continuous or because it contains frequencies other than those in the 

1 5 basis of the output of the Fourier transform, significant errors are produced. These errors 

manifest themselves as a spreading of the spectrum in which amplitude (or po\yer) is spread out 
over frequencies that are not actually in the signal. Fourier methods approximate the actual 
frequency content of the signal by spreading amplitudes to frequencies away from the basis 
frequencies. Such spreading is a well-known problem with Fourier methods and is called 

20 spectral broadening. Spectral broadening cannot be eliminated, although various methods have 
been employed to reduce it such as those that rely on windows. These windowing methods are 
not relevant to the current invention other than their well known existence and extensive use 
confirms that spectral broadening is a significant problem caused by the discrete approximations 
that are inherent in Fourier methods. Spectral broadening will always occur when discrete 

25 frequencies are used to estimate the spectral content of signals that have spectra that are either 
continuous or composed of such a inultitude of closely spaced frequencies that they are 
essentially continuous, such as the signals generated by Doppler shifting of an interrogating 
signal from blood flowing in an artery. 
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Another major problem of Fourier methods is that they cannot be used to examine all 
frequencies of interest. They can only include frequencies that are harmonics of the fundamental 
frequency used in the measurement. These problems, again, are well known and they limit the 
ability of Fourier methods to examine frequencies that are lower than the fundamental frequency 
5 used in the measurement. Furthermore, this limitation leads to the frequency resolution equaling 
the fundamental frequency. For example, to estimate the amplitudes of frequencies spaced 20 Hz 
apart the fundamental frequency must be 20 Hz or lower. Therefore, if an interrogating signal 
has a frequency of 2 MHz and the need exists to discriminate the received signal with a 
resolution of 20 Hz than at least 100,000 harmonics are in the basis. This large number of 
10 harmonics creates significant computational difficulties and is generally regarded as impractical 
for a variety of reasons, including the requirement to use many data points to perform such 
calculations. 

The limitation that all basis frequencies for Fourier series are based on a fundamental 
frequency also causes other problems, such as when the received signal is not stationary (varies 

15 with an independent variable, such as time). This situation will exist, for example, when the 
ireceived signal has its spectral content changing as the result of the interrogating signal 
interacting with blood flowing through an artery. The velocity profile changes in an artery as the 
blood pressure varies in response to the beating of the heart and the variations in the blood's 
velocity profile will lead to a nonstationary signal that has a time-varying SD function. If the SD 

20 function is to be measured over a short time interval, such as over 0.01 second, then the 

fundamental frequency for the Fourier analysis cannot provide a better frequency resolution than 
that which is associated with the measurement time interval. For example, a data sample with a 
period of 0.01 second will have a fundamental frequency of 100 Hz, which precludes calculating 
the spectral content of signals with a frequency resolution better than 100 Hz . 

25 Paranietric methods also suffer deficiencies that are significant for certain applications. 

In particular, these methods minimize the sum of the squared errors between the measured data 
and the values estimated by the model. As described later, measured data are necessarily 
digitized and this process introduces errors that produce a wide variety of equally valid models. 
Minimizing least squared errors does not account for the possibility of such alternative solutions. 
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Another shortcoming of parametric methods is their critical dependence on the order of 
the model chosen. For example, in the most common model, AR models, if too few poles are 
used then spectral peaks are not reliably detected. If too many poles are selected then false peaks 
are introduced. Proper model selection is nontrivial and requires a computationally intensive 
5 search process that involves calculating models using a range of orders and then selecting among 
them. 

It has also been recognized that digitization errors can significantly affect spectral 
analyses in certain applications. The input analog signal is not used directly to calculate spectral 
density functions. First, the measured data are digitized and given discrete values. This process 

10 quantizes the data and introduces errors. These digitization errors (also called discretization 
errors and quantization errors) cohfound spectral analysis. Each digitized measurement is an 
estimate of the actual value and has an upper and lower error limit, within which any value is 
equally likely to occur. Even though any value within the error limits is equally likely to be the 
true value, the digitized value is conventionally established as the midpoint of the range between 

15 the lower and upper error limits. Any value within the range is equally likely, leading to an 

infinite number of possible values that can be used to calculate spectra. Thus an infinite number 
of candidate spectral density functions exist. Existing spectral analysis methods do not 

explicitly account for this range of possible solutions. These methods can be interpreted as using 
least squared errors or other niinimum norm methods that attempt to minimize differences 

20 between predicted values and measured values. The structure of these problems makes them 
sensitive to such errors when they are solved using minimum norm methods, and in particular 
when solved using least squared errors methods and produce an ill-posed problem. These ill 
posed problems are well known and regularization techniques are used to calculate reasonable 
estimated solutions. Regularization methods, including Tikhonov regularization (sometimes 

25 called ridge regression or Phillips regularization) are used with least squared errors fits and 
augment the least squared errors model with an auxiliary function that is extremized (such as 
minimized). Regularization methods produce stable results, but trade off quality of the fitted 
solution (comparison of the predicted values to the measured data) with some measure of the 
regularization function. Making this tradeoff requires using a weighting factor, which is also 



called the regularization parameter, that trades off the weight given to fitting the data with the 
weight given to the regularization function. 

Existing methods for solving ill-posed problems have several deficiencies. They require 
calculating the weighting factor, which is difficult and can have a significant effect on the results. 
5 They fail to explicitly recognize that a solution is only valid when all of the estimated values are 
within the digitization error bins of the measured values. This deficiency manifests itself in 
solutions that have estimated values that are outside of the digitization bin for one or more 
measured values, and, thus, are not correct. Another deficiency is that in practice the methods 
rely on using minimizing least squared errors, which tends to force estimated values toward the 

10 center of the digitization bins, thus not providing a monotonically convergent means for 
calculating viable solutions by adjusting the regularization parameter. 

Therefore, significant deficiencies exist in the prior art regarding the sjpectral analysis of 
time-varying signals and signals when it is desired to obtain fine frequency resolution. 

The present invention overcomes problems associated with determining spectral content 

15 from signals where the fundamental frequency of the signal being processed corresponds to a 

time period longer than what has been measured and/or the underlying spectrum is, at least over a 
range of interest, continuous or substantially continuous (occasionally hereinafter collectively 
referenced as "substantially continuous")- The present invention fixrther enables spectral analysis 
where the signal is nonstationary and where the signal has digitization errors. 

20 More specifically, the present invention is directed to a method and apparatus ("utility") 

involving a time-to-frequency domain transform to define a substantially continuous spectrum or 
other spectrum including non-zero values at irregularly spaced frequency intervals. The invention 
thus allows for spectral estimates for more than a small number of discrete frequencies thereby 
enabling a more accurate representation of certain real spectra and reducing or eliminating 

25 spectral broadening. The invention also allows for examination of a greater range of frequencies 
of interest and provides frequency resolutions that are not limited to a fundamental frequency. In 
this regard, the invention allows for determination of a spectrum based on analysis of a portion of 
the time-based signal that is shorter than a full period of a component of the time-based signal. 
Moreover, the present invention allows for enhanced analysis of received signals that are not 
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stationary, i.e., that vary with an independent variable such as time. Therefore, the invention has 
particular advantages with respect to medical contexts where the received signal has a changing 
spectral content as a result of interaction of an interrogating signal with moving physiological 
material such as blood flowing through an artery. The invention also enables reduction in 
digitization or quantization errors associated with analog-to-digital conversion. 

According to one aspect of the present invention, a utility is provided for processing an 
input signal, such as a medical diagnostic signal, to obtain a spectrum that includes non-zero 
values at irregularly spaced intervals. An associated process involves receiving time-based 
information corresponding to one or more predetermined time intervals of a time-based signal, 
performing a transform on the time-based information to obtain a frequency spectrum including 
non-zero values at irregularly spaced intervals, and operating a processor in a signal processing 
environment for using the transform to provide an output based on the time-based signal. The 
spectrum includes at least one pair of successive points (i.e. non-zero values) that has a different 
frequency spacing than at least one other set of successive points. Preferably, the spectrum 
includes a set of non-zero values including a first nonzero value at a first frequency and a second 
non-zero value at a second frequency greater than the first frequency where the second fi-equency 
is a non-integer multiple of each fi-equency of the set of frequencies other than the second 
frequency. 

According to another aspect of the present invention, a utility is provided for processing 
an input signal, such as a medical diagnostic signal, so as to generate a spectrum based on 
information corresponding to a short time interval of the input signal. A corresponding method 
includes the steps of receiving time-based information corresponding to one or more defined time 
intervals of the time-based signal, where the time-based signal includes a component having a 
period that is longer than the time interval, performing a transform on the time-based information 
to obtain a frequency spectrum for the time-based signal, and operating a signal processor in a 
signal processing environment for using the transform to provide an output based on a time- 
based signal. This process can be used to obtain a frequency spectrum where the received time- 
based information corresponds to a time interval that is less than half of the period of the signal 
component. 
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According to a further aspect of the present invention, a utility is provided for processing 
an input signal to obtain a substantially continuous spectrum. An associated method includes the 
steps of receiving time-based information corresponding to one or more defined time intervals of 
a time-based signal, performing a transform on the time-based information to obtain a frequency 
5 spectrum for the time-based signal where the spectrum defines a substantially continuous 

function across a fi'equency range that has non-zero values for a majority of firequencies in the 
range, and operating a processor in a signal processing environment for using the transform to 
provide an output based on the time-based signal. Such a substantially continuous spectrum can 
more accurately characterize the real spectrum under analysis, thereby avoiding spectral 

1 0 broadening and other shortcomings of certain prior art processes. 

According to a still further aspect of the present invention, a utility is provided for 
processing an input signal so as to account for a digitization or quantization error associated with 
analog to digital conversion. An associated method includes the steps of receiving time-based 
information corresponding to one or more defined time intervals of a time-based signal, where 

1 5 the time-based signal is an analog signal and the time-based information is digital information, 
performing a transform on the time-based information to obtain a frequency spectrum for the 
time-based signal, where the step of performing a transform involves accounting for a 
digitization error associated with a difference between the analog time-based signal and the 
digital time-based information, and operating a processor in a signal processing environment for 

20 using the transform to provide an output based on the time-based signal. The invention thus 
allows for accurate definition of a spectrum based on the input signal for certain cases where 
accurately defining such a spectrum may be difficult or impossible due to digitization error. 

In accordance with another aspect of the present invention, a utility is provided for use in 
connection with converting a signal from an analog form to a digital form. The utility involves 

25 receiving an analog signal, defining a number of digitization ranges or bins, and assigning a 
specific value within a bin to a sample related to the analog signal based on a mathematical 
model where the specific value can vary relative to the boundaries of the bin. Preferably, the 
mathematical model determines the specific value by solving a constrained optimization problem 
involving an objective function having one or more constraints. For example, such constraints 
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may relate to one or more limit values of a bin, a nonnegativity constraint, and a peak count 
within the bin. For example, the objective function may be constrained to be convex and positive 
within the bin range with slopes determined by predefined parameters at or near the bin limits. 
The constraints may be implemented as penalty or barrier functions, for example, Heaviside 
functions. In one implementation such a function approximates a Heaviside function and is 
tapered at an area corresponding to a constraint value (e.g., a bin limit value) such that the 
function is free from singularities (undefined function values) at that area. 

In one implementation, the invention is used to determine dimension-related information 
regarding a flow channel of a physiological fluid based on analysis of a signal modulated based 
on interaction with the fluid of the flow stream. A corresponding system includes a detector for 
receiving a signal from the fluid and providing an analog electronic signal based on the received 
signal, an analog-to-digital convertor for providing a digital signal based on the analog electronic 
signal, a processing means (processor) for first processing the digital signal to obtain spectrum 
information corresponding to a frequency spectrum including non-zero values that are not at 
regularly spaced frequencies, and a second processing of the spectrum information to obtain 
dimension-related information for the flow channel. The dimension-related information may be, 
for example, a channel dimension or a volumetric flow rate of the fluid. 

BRIEF DESCRIPTION OF THE FIGURES 

For a more complete understanding of the present invention and further advantages 
thereof, reference is now made to the following Detailed Description taken in conjunction with 
the figures as described below: 

FIG. 1 Example Spectral Density Function 

FIG. 2 Illustrates an example Gaussian pulse 

FIG. 3 Example Amplitude Function Made from 10 Piecewise Continuous 

Functions 

FIG. 4 Example Amplitude Function Composed of Multiple Functions with 

Overlapping Domains 
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FIG. 5 Example Phase Angle with Doppler 

FIG. 6 Example Phase Angle with Doppler when Emitted Frequency Selected to 

Produce Approximately Linear Phase Angle Response 

FIG. 7 Example Phase Angle with Doppler when Emitted Frequency Selected to 

5 Produce Always Positive and Approximately Linear Phase Angle 

Response 

FIG. 8 Example sin(Phase Angle) with Doppler when Emitted Frequency 

Selected to Produce Positive Phase Angle Response Close to Zero Phase 
Angle 

1 0 FIG. 9 Example cos(Phase Angle) with Doppler when Emitted Frequency 

Selected to Produce Positive Phase Angle Response Close to Zero Phase 
Angle 

; FIG. 1 0 Example Phase Angle with Doppler Using 2 MHz Emitted Signal 

FIG. 1 1 Example cos(phase angle) with Doppler Using 2 MHz Emitted Signal 

15 FIG. 12 Example sin(phase angle) with Doppler with 2 MHz Emitted Signal 

FIG. 13 Six Parameter Distribution for Velocity 

FIG. 14 Quadratic Functidn Segment Expression for f(t) without Continuity 
FIG. 15 Parameters for Example 
FIG. 16 Example ||A|| Matrix 
20 FIG. 17 Example ||P|| Matrix 

FIG. 18 Example Envelope Function 

FIG. 19 Spectral Density Component Functions: Used to Calculate Example f(t) 
values 

FIG. 20 Example Signal and Locations of Measurement Points 
25 FIG. 21 Example Signal after Adjust for Gaussian Envelope and Locations of 

Measurement Points 
FIG. 22 Exariiple Calculated Results for G(f) and Actual G(f) 
FIG. 23 Example Calculated Results for F(f) and Actual F(f) 
FIG. 24 Equations for piecewise Continuous Equally Spaced Quadratic Function 
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Segments 

Equations for Piecewise Continuous Equally Space Linear Function 
Segments 

Three Frequency Function with Four Digitization Bins 

Three Frequency Function with Four Digitization Bins and Least Squared 

Errors Estimate 

Three Frequency Function with Four Digitization Bins and Least Absolute 
Value Errors Estimate 

Three Frequency Function with Four Digitization Bins Least Squared 
Errors and Least Absolute Value Errors Estimates 
Actual and Estimated G(f) Calculated without Accoxmting for Digitization 
Actual and Estimated F(f) Calculated without Accounting for Digitization 
Bounded Area Measure when Using Piecewise Continuous Quadratic 
Function Segments 

Arc Length Measures for Spectral Density Component Functions Based on 

Piecewise Continuous Quadratic Function Segments 

Bounded Area Measures for Spectral Density Component Functions Based 

on Piecewise Continuous Quadratic Function Segments 

Quadratic Curvature Measures for Spectral Density Component Functions 

Based on Piecewise Continuous Quadratic Function Segments 

Objective Function Combining Measure of Arc Length and Measure of 

Area for both G(f) and F(f) 

Double Sided Constraint Violation Value Function 

Nonnegativity Constraint Value Function 1 (NCVl) 

Nonnegativity Constraint Value Function 2 (NCV2) 

Tapered Heaviside Gate Function Examples 

Tapered Heaviside Gate Function for Extreme Value of G(f) Less than 
Zero 

Example G(f) Calculated Using NCVl Nonnegativity Constraint Value 
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Function 

FIG. 43 Peak Count Constraint Value (PCV) Function 
FIG. 44 Basic Calculation Sequence (Page 1 of 3) 
FIG. 45 Basic Calculation Sequence (Page 2 of 3) 
FIG. 46 Basic Calculation Sequence (Page 3 of 3) 
FIG. 47 Actual and Calculated G(f) from Digitized Data (M=50) 
FIG. 48 Actual and Calculated F(f) from Digitized Data (M=50) 
FIG. 49 Actual and Calculated Spectral Density Function from Digitized Data 
(M=50) 

FIG. 50 Actual and Calculated G(f) from Digitized Data (M=200) 
FIG. 5 1 Actual and Calculated F(f) from Digitized Data (M=200) 
FIG. 52 Actual and Calculated Spectral Density Function from Digitized Data 
(M=200) 

FIG. 53 Actual and Calculated Cumulative Spectral Density Function from 

Digitized Data (M=200) 
FIG. 54 Block diagram of one embodiment of a system in accordance with the 

present invention 

DETAILED DESCRIPTION 

The interrogating signal emitted from an energy emitter, such as a transducer, can take on 
many forms and EQ. 1 is a general representation of an interrogating signal with frequency / 
f(0 = E(Osin(2 7C/0 (1) 

E(t) is a fiinction that describes how the amplitude varies with time. For example, if the 
emitted signal is a pulse then E(t) can be a fiinction that includes a phase shift term that will 
periodically repeat with period x. An example of a signal that has such a form is a Gaussian 
pulse. The amplitude of a Gaussian pulse will vary with time as defined by EQ. 2 
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(2) 
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A Gaussian pulse will vary with time as defined by EQ. 3 and a single Gaussian pulse is 
illustrated in FIG. 2. 



( 2 2^ 



.2 

, 1 4^e^ UmilKft) (3) 



These equations can be used with either pulsed or continuous wave systems. In the 
context of ultrasound medical analysis, the pulsed wave has the advantage of allowing selection 
of the tissue depth to interrogate, which can reduce errors caused by Doppler shifts caused by 
blood deeper than the region of interest, such as the ascending aorta. The pulsed wave approach 
leads to pulses that have a transient response, and this needs to be dealt with, possibly by waiting 
for the transient to pass out of range before collecting data. Another approach could be to model 
the transient explicitly as a pure tone that has changing amplitude or as a pure tone that has 
multiple amplitudes that phase shift with time. Using a continuous wave avoids the transient 
response issues and leads to a pure tone being used to interrogate the tissue. The continuous 
wave approach, however, has the confounding factor of receiving Doppler-shifted signals from 
regions other than the selected depth of the ascending aorta, or whatever other region is of 
particular interest. Depth dependent signals could be inferred from signals that are created from 
different emitted frequencies (these would need to be much different, such as 1 MHz and 3 MHz) 
such that attenuation effects are substantially different. The signal strength differences could be 
matched against spectra to start getting depth-related data. Higher frequencies attenuate more 
rapidly than lower ones so using higher frequencies reduces the interference from deeper tissues, 
but at the expense of reducing the signal from the region of interest. 

As described earlier, the interrogating signal can be used to energize tissue. The resulting 
backscattered signal can be sensed and analyzed to determine useful properties of the tissue that 
backscattered the signal. Such a signal can comprise, for example, the received signal that echos 
back from a multitude of reflectors or backscatterers moving at one or more velocities that have 
been insonified by an ultrasonic interrogating signal. Among the purposes for which the present 
invention can be used is determining the SD functions of backscattered signals from interrogating 
signals that have either constant or time-varying amplitudes, including time varying signals that 
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include one or more pulses, repeating sets of pulses, or any other amplitude function that varies 
with an independent variable. As discussed above, time is an example of such an independent 
variable. 

For illustrative purposes, the following description and derivations will use as an example 
signals such as those that may occur when ultrasonic signals are used to interrogate tissues in a 
living organism. This illustration is not limiting in that the present invention does not require 
using interrogating signals or a response from such a signal. The measured data used by the 
present invention may be acquired using any suitable means that eventually produces two or 
more digitized data, including transducers followed by subsequent processing or data input or 
read from one or more sources of data that may have been previously stored, such as tables or 
graphs stored as hard copy or using electronic, magnetic, or optical means. 

Let A(t J) be the SD function (in this case the function is the same as an ampHtude density 
fimction). The general form of the signal for which the SD function is to be obtained is EQ. 4. 



As described in more detail below, the present invention determines SD functions for 
processes with two or more independent variables, such as signals that vary with time, g(t), when 
the SD function depends on at least one independent variable,/, that has one or more nonzero 
values over one or more parts of the interval [fi^fn ]j where and // are the lower and upper 
limits of f. EQ. 4 defines an amplitude function, A(t,J), that varies based on the values of two 
independent variables, t and / For this example, the variable t represents time and the variable / 
represents frequency. The A(t, f) function can be separated into functions of f and of t. 
A(/,/) = E(/)B(/) (5) 

When EQ. 5 is substituted into EQ. 4 the result is a general expression for the time 
varying signal g(t) for which the SD function B(f) is sought. 

g(0= f E(OBC/)sin(2 7i/r + (])(/))# (6) 



H 
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EQ. 6 can be rearranged to become EQ. 7 

■ I" 

g(0= E(OB(/)(sin(2K/Ocos((|)(/)) + cos(2 7i/Osin(<|)(/)))t// (7) 

\ 

The following terms are defined for coefGcients. 
F(/)=B(./)sm(<|)(/)) (8) 
G(/)=B(/)cos((t)(/)) (9) 

Substituting these equations for the coefficients into EQ. 7 leads to EQ. 10 

I" 

g(0 = E(0 F(y)cos(2 7c/0 + G(/)sin(2ji/0^// (10) 
\ 

Distributing the integral over the sum of the sine and cosine terms leads to EQ. 1 1 . 




Define spectral density component functions (SDC functions) as functions from which 
the SD function can be calculated. F(f) and G(f) are SDC functions. As is well known, F(f) and 
G(j) can be used to calculate the spectral density function B(f) using EQ. 12. 

B(/)= Vf(/)' + G'(/)^ (12) 

The amplitude functions can also be used to calculate the phase shift for each frequency 
using EQ. 13. 

(^(y) = arctanj^|^j (13) 

SD functions and SDC functions both characterize how the amplitude of a function varies 
with respect to an independent variable and are collectively referred to as amplitude functions. 
For example, F(f) characterizes how the amplitude of F(f) cos(27rf t) varies with frequency. 
Either the SD function itself or SDC functions are to be estimated using the means described 
below. 

The SD function or the SDC functions, such as F(f) and G0), may be represented using 
various functions. The functions used to represent amplitude functions will generally have 
parameters that, when selected correctly, will lead to useful estimated values that are calculated 
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for g(t). Amplitude functions may be represented by any functional form for which parameters 
can be determined that produce estimated values for g(t) that differ from the actual values of g(t) 
in a manner that meets predetermined objectives. An example of such a predetermined objective 
is having the errors between the measured values and estimated values of g(t) being small enough 
that the estimated values for g(t) provide useful approximations of the actual values of g(t). 
Useful values for g(t) occur when the parameters specified for the functions that represent the 
SDC functions can be used to calculate the SD function of g(t). 

One aspect of the present invention involves structure and methodology for using values 
of g(t) and E(t), either based on measurements, simulations, or other suitable means, to determine 
the specific values of the parameters that characterize the SD function or one or more SDC 
functions when such amplitude functions are composed of one or more substantially continuous 
functions or when such amplitude functions use values for one or more independent variables, 
such as frequency, that lead to non-orthogonal basis functions constituting g(t). Using 
continuous functions to represent one or more amplitude functions over at least part of the 
domain of an independent variable (such as time) is particularly advantageous when a multitude 
of signails with various frequencies combine to produce a composite signal that is the measxu'ed 
signal g(t). Such a composite signal occxurs, for example, when an interrogating signal 
backscatters from the multitude of particles traveling over a wide range of velocities as they are 
being carried in a flowing fluid. An example is the signal produced when an interrogating 
ultrasonic signal backscatters from the blood cells traveling through a living organism's aorta. 

hi contrast to the present invention, Fourier methods represent the amplitude functions 
using discrete values of frequencies that aire selected as integer multiples of a fundamental 
frequency. The frequencies used in Fourier methods are carefully selected such that the resulting 
basis functions are both orthogonal and normal (orthonormal). As described earlier, using 
discrete frequencies, particularly when they are restricted to being harmonics of a fimdamental 
frequency, to approximate continuous functions leads to significant inaccuracies, such as spectral 
broadening. The present invention improves upon Fourier methods by using frequencies that can 
be selected substantially without limitation as to which frequencies are selected and also by 
allowing the use of continuous functions to represent continuous amplitude functions. 
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In the present invention, one or more functions are selected to represent the SD or SDC 
functions. The function used to represent an SDC, such as F(f), may be zero. The functions used 
to represent the ampHtude functions may be discrete, such as one or more vectors of discrete 
frequencies starting at any frequency of interest and spaced at any, possibly varying and 
nonuniform, frequencies from each other (in contrast to the evenly spaced harmonic frequencies 
using in Fourier methods). The functions representing amplitude functions in the present 
invention can be orthogonal or nbn-orthogonal. The functions representing amplitude functions 
can be either normal or not normal. The function or functions used to represent amplitude 
functions may be one or more continuous, not continuous, or piecewise continuous functions. 
The functions representing amplitude functions may have one or more intersecting, including 
overlapping, domains. An overlapping domain between two or more fiinctions is an interval 
where the domains of two or more functions intersect at more than one value. FIG. 3 illustrates 
an example where multiple piecewise continuous functions that intersect at only their endpoints 
represent an amplitude function. FIG. 4 illustrates an example where multiple functions with 
overlapping domains model an amplitude function. One aspect of the present invention is that 
functions representing the amplitude functions may use one or more of the following functions: 
splines, B-splines, polynomials of order 1 or greater, trigonometric functions, trigonometric 
polynomials, exponential functions, hyperbolic functions, Heaviside functions, functions that 
approximate Heaviside functions (discussed in detail below), or any another functions that can be 
used to construct one, or more bases that span one or more predetermined vector spaces that can 
approximate or define a discrete, non-continuous, continuous, or piecewise continuous function. 
For example, the piecewise continuous segments illustrated in FIG. 3 are composed of quadratic 
polynomials and Heaviside functions. As another example, FIG. 4 illustrates how B-splines can 
combine to represent an amplitude function. 

FIG. 4 also illustrates another aspect of the present invention that provides for 
computational flexibility. The independent and dependent variables can be offset (as illustrated 
in FIG. 4 where the frequency axis is offset by 2 MHz) or scaled to facilitate calculating the 
parameters of the functions that represent amplitude functions. 

Using piecewise continuous functions or functions with overlapping domains is of 
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particular interest because parameters for function segments defined over discrete intervals can 
be efficiently determined using the present invention. The union of the function segments can 
represent complicated amplitude functions. FIG. 3 is a graphical example of a piecewise 
continuous function. The points where the . individual function segments of the piecewise 
continuous function intersect are knot points and these are illustrated by the hollow dots in 
FIG. 3. 

Another aspect of the present invention involves structure and methodology employing a 
calculation module that uses one or more objective measures to determine the values of the 
parameters in the functions that represent the amplitude functions or SD function. One example 
of such an objective measure is using goodness of fit. Other examples of objective measures are 
minimizing the length of the amplitude functions, SD function, or SDC functions, minimizing 
the squared area or absolute value of the area under the amplitude functions, SD function, or 
SDC functions, and minimizing one or more measures of the slopes or changes in slopes of the 
amplitude functions, SD function, or SDC functions. 

One example of an objective measure using goodness of fit is minimizing, or 
approximately minimizing, the sum of the squared errors between measured values of g(t) and 
those predicted using amplitude functions. Another example of an objective measure of 
goodness of fit is minimizing, or approximately minimizing, the sum of squared errors when the 
values of one or more of the parameters in one or more of the amplitude functions or the SD 
function are constrained to have values with predetermined properties. Such predetermined 
properties may include that the parameters are allowed to only have particular values, such as 
being non-negative, or that the parameters have a particular relationship, such as the 
mathematical combination of two or more of the parameters (such as one or more of their sum, 
difference, product, or quotient,) always be positive or that the mathematical combination always 
falls within a certain range, such as being approximately equal to zero. Another exatnple of an 
objective measure of goodness of fit is minimizing, or approximately minimizing, the sum of the 
absolute values of errors, including the possibility that the values of one or more of the 
parameters in one or more of the amplitude functions or the SD function are constrained to have 
values with predetermined properties. As before, such properties include having one or more 
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parameters being restricted to particular values, such as being non-negative, or that the 
parameters have a particular relationship, such as the sum of two or more parameters always be 
non-negative or that they approximately equal zero. These aspects of the invention are covered 
in more detail below. 

5 Referring to EQ. 1 1 , in accordance with a further aspect of the present invention, 

optimized results can beobtained by setting the values for^ and respectively low enough and 
high enough that they cover the full range of frequencies that constitute the spectrum of g(t). The 
values for fi^ and may be determined concurrently with determining the values of the 
parameters for the functions representing the SDC functions, such as G(/) and F(/), or they may 

10 be selected prior to determining the parameters for the functions representing G0 and F(/). 

Selecting^ and ahead of time is often possible. For example, the conditions that produce g(t) 
are generally known. A predetermined lowest frequency may be based on the frequency of ah 
interrogating signal interacting with one or more moving media. For example, the lowest 
frequency can be predetermined based on the frequency of the interrogating signal and the lowest 

15 frequency expected from the moving media after including factors such as measurement errors or 
noise that can affect the frequencies produced by the interaction of the interrogating signals with 
the moving media. 

Another aspect of the invention involves using single continuous functions that span the 
entire domain of frequencies as amplitude functions or as SD functions. An example of such a 
20 function is EQ. 14. This equation characterizes the amplitude functions that occur during 
laminar fluid flow with measurement noise. This equation would become B0 in EQ. 8 and 
EQ. 9. 
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The emitted frequency sent out as the interrogating signal is and the tissues move at 
various velocities to produce the reflected frequencies,/ The amount of tissue at each velocity 
produces the relative amplitudes. The parameters are: 

1. F the fraction of the insonified region that has flowing liquid. 

2. k: the shape exponent for the flow region 

3. : the standard deviation for the velocity in the slow/nonflowing region's reflectors. 

4. : the frequency of the interrogating signal. 

5. a: the maximum value for velocity that is' achieved at the time that a measurement is 

made. 

6. ^i^ : the mean velocity of the slow moving/nonflowing region's reflectors. 

7. : the standard deviation of the measurement error associated with the flowing fluid. 

8. c: the speed of sound. 

A further aspect of the present invention involves using such a single continuous 
amplitude function in combination with a function that represents the phase shift, ip(/), see, e.g., 
EQ. 8 and EQ. 9, to produce single continuous functions that represent amplitude functions that 
span the entire domain of frequencies. An additional aspect of the present invention involves 
using such a single continuous amplitude function, e.g., in combination with functions that 
represent the sin(4>(f)) and cos(4>(f)) in EQ. 8 and EQ. 9 to produce single SDC functions that 
span the entire domain of frequencies. A fxirther aspect of the present invention relates to 
representing ^(f), sin(<f>(f)), or cos((p(f)) using either constants or polynomials with degree of 1 or 
greater. A further aspect of the present invention involves selecting the frequency of the 
interrogating signal used with Doppler measurements of velocities such that <l>(f), sin(4>(f)), or 
cos(<t>(f)) may be represented using either constants or a linear function. 

A further aspect of the present invention involves characterizing the phase shift when 
using Doppler measurements of fluid flows using EQ. 15 to represent phase angle. 
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where 

c|) = phase angle in radians 
f = frequency in Hz 

D = depth of region from which Doppler-shifted signals are being generated 
c = speed, of signal, such as speed of sound 
X = attenuation exponent 
b = attenuation constant 

For the case where x=l EQ. 15 simplifies to become EQ. 16 
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The attenuation constant, b, can be calculated from the attenuation factor using EQ. 17 
b ^ .5000000000 10 p ln( 10 ) (17) 

where P is the attenuation factor in dB. 

FIG. 5 illustrates how the phase angle can change with frequency. As is known, selecting 
the emitted frequency determines the frequencies of the returned signal. As previously 
mentioned, one aspect of the present invention involves selecting the emitted frequency such that 
the returned signals are confined to an approximately linear region of the phase arigle response. 
For example, in FIG. 5 selecting an emitted firequency that leads to retum fi-equencies around 
2.079 MHz would be in a region where phase angle varies close to linearly with frequency, 
whereas selecting an emitted fi-equency that leads to return frequencies about 2.17 MHz would be 
in a region where phase angle does not vary linearly with frequency. More specifically, if 2.078 
MHz is selected as the emitted frequency and blood flow in the ascending aorta of an adult 
human is being measured using an interrogating signal beamed down the axis of the aorta, then 
firequencies in the range of approximately 2.078 to 2.093 MHz will be returned. FIG. 6 illustrates 
an example phase angle response over this fi-equency range. In this example, selecting 2.078 
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MHz leads to phase angle responses that may be modeled as a linear function. 

FIG. 5 illustrates that half of the phase angle response will be over frequencies that lead 
to positive phase angles. Selecting an emitted frequency such that the phase angle response is 
always positive is beneficial. When the phase angles are always positive the amplitude response 
function of EQ. 1 will always be positive. If the phase angle becomes negative the sin(4>(f)) 
term will be negative, leading to negative values over at least part of the range ofF(f). Having 
only positive phase angles avoids this condition. Calculating correct values for F0) and G(f) is 
made easier when it is known that F(f) and G(f) are always positive. Selecting the emitted 
frequency for Doppler velocity measurements such that the phase angle response is always 
positive is included in accordance with another aspect of the present invention. 

Another aspect of the present invention involves selecting the emitted frequency so that 
the phase angle response is always positive and in a region where the phase response is 
approximately linear. FIG. 7 illustrates an example in accordance with this aspect of the present 
invention. 

A further aspect of the present invention involves selection of the emitted frequency so 
that at least a portion of the phase angle response is close to zero. 

FIG. 8 and FIG. 9 illustrate how thus selecting the emitted fi'equency yields sin((p(f)) and 
cos(<t>(f)) functions that can be approxiniated. sin(4>(f)) can be approximated in many cases using 
a linear function. The range of cos((f>(f)) is small. Depending upon the precision needed it may 
be approximated as a constant (such as 0.98 in this example), as a linear function, or as a higher 
order polynomial. A linear approximation would pick up most of the variation and be adequate 
for many applications. 

Compare the results illustrated in FIGS. 7, 8 and 9 with what happens if the emitted 
fi'equency is not selected according to this aspect of the present invention. FIGS. 10, 11, and 12 
illustrate the results from using an emitted frequency of 2 MHz for the conditions of this 
example. The linear response and small domains that occurred when the emitted frequency was 
selected in accordance with this aspect of the present patent are not present in FIGS. 10, 11, and 
12. 
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When the frequency distribution of a spectral distribution has been produced in a manner 
that relates frequency to a calculated variable then the distribution of the calculated variable can 
be determined. For example, if the frequency distribution has been obtained using Doppler then 
the velocity at any frequency can be calculated using EQ. 18. The frequency data or data for a 
calculated variable can be used to calculate the values of one or more parameters that 
characterize a distribution. Such distributions are particularly useful when the methods of the 
present invention are used to calculate continuous distributions. Such a distribution could be, for 
example a Gaussian distribution that is characterized by its mean and standard deviation. EQ. 19 
is an example of a more complicated distribution that has six parameters. FIG. 18 is an example 
of such a distribution when the calculated parameter is velocity. 
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Other aspects of the present invention involve structure and methodology for determining 
the parameters of an SD fiinction, or the SDC fimctions, from signals that have been distorted by 
digitization errors. When signals are measured and converted from analog to digital form the 
analog values are converted to numeric (digital) values with a limited precision. Such limited 
precision leads to numeric values that approximate the actual values of the signal and differences 
between the actual values and the approximate values are errors caused by the digitization 
process. In particular, the present invention includes structure and methodology for calculating 
SD fimctions, or the fiinctions that can be used to determine the SD fimctions, when using 
digitized data by constraining the calculated fiinctions so that they have predetermined properties. 
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Such predetermined properties can include having at least one, and typically the majority or 
approximately all, of the values of g(t) estimated using the calculated functions falling within the 
ranges determined by the digitization process. Such predetermined properties can also include 
requiring SD functions, or the functions that can be used to determine the SD functions, having 
predetermined properties related to their geometries or of the parameters that define the 
functions. For example, such methodology includes having functions that are the minimum 
measures of slopes or rates of changes in slopes, the lengths, or squared areas enclosed by the 
functions. Another example is having functions that maximize, or minimize, the magnitudes of 
one or more of the parameters that characterize the functions. Another example is having 
functions that minimize measures of length and squared areas enclosed by functions and that 
maximize the magnitudes of one or more of the parameters that characterize the functions. These 
aspects of the invention are described more completely, later. 

To illustrate the aspect of the present invention that characterizes amplitude functions 
using multiple functions, consider the following description of a general process for estimating 
amplitude functions. Amplitude functions, such as those of EQ. 8 and EQ. 9, can be 
approximated arbitrarily closely as the sum of other functions. The following expressions define 
the approximations of G(f) and F(f) as the sums of ATpiecewise functions. 

0{f)=J,Gff) (19) 

F(y)=Z^;(/) (20) 

7=1 

The ATpiecewise functions defined over each of N intervals will collectively span the 
domain of f that is of interest. N functions exist for G(f) and N functions exist for F(f), and they 
may or may not have overlapping domains. Each of the Gj(f) and Fj(f) will be defined over a 
domain of/ with lower and upper limits. For the case of using piecewise continuous functions 
the union of the domains will be continuous arid the individual domains will not overlap except 
at the points of intersection. 
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Substituting EQs. 19 and 20 into EQ. 1 1 and distributing cosine and sine terms over the 
sums leads to EQ. 2 1 . 
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Changing the integral of the sum to a sum of integrals leads to EQ. 22. 
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Each of the piecewise functions is defined over part of the domain of/ Indexing the 
limits of integration for each of the integrals in the sums leads to EQ. 23. 
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E(t) is the envelope of the interrogating signal. If the signal being analyzed has not been 
modulated by an interrogating function then setting the envelope function E(t)=l will decompose 
the received signal without removing the effects of the interrogating signal. If E(t) is not constant 
then methodology in accordance with the present invention can separate the effects of E(t) while 
decomposing the received signal into, for example, SD or SDC fimctions. 
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g(t) is the received signal and measured values for it are known. Dividing both sides of 
EQ. 23 by E(t) produces EQ. 24. 
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Measured values from g(t) are the data used to calculate F(f) and G(f). 
Use the following substitutions in EQ. 24.: 

g(0 



E(/) 



= f(0 (25) 



^H.-^j (26) 

4=^y (27) 
and obtain EQ. 28. 
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As stated earlier, the piecewise continuous fimctions may be polynomials, such as linear 
or quadratic functions. Each of the function segments will have its coefficients, such as the slope 
and intercept for linear functions. For illustration of the principles of the present invention, the 
following example uses quadratic functions and derives expressions that may be used to solve for 
the coefficients that specify each of the quadratic function segments. 

Let the following expressions specify the piecewise quadratic function segments: 
Gff)^a/^h,f^c, (29) 
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P'P)=x.f+y/+Zj (30) 

Substituting these expressions into EQ. 28 leads to EQ. 31. 



f(0 = 
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This equation is linear in the coefficients Xp Zy, and Cy, and the integrals have 
analytical solutions. FIG. 14 has the analytical solution. This equation is a general expression. 
No restrictions have been placed on the interval limits [Ly, so they do not need to be 
continuous and they could overlap. Therefore, the function segments could be discontinuous. No 
constraint has been placed on the slopes of function segments matching if function segments are 
continuous. The current formulation, without constraints, is broadly applicable. 

The unknowns to be solved for are the coefficients jCy, y^^ Zp Up bp and Cj for j-1 . . N. 
Represent these unknowns as the vector |d|. Once values have been estimated for |d| then 
estimated values for f(t) can be calculated using EQ. 3 1 . The values of |d| can be solved for using 
one or more objective criteria, such as the estimated values for f(t) having the least sum of 
squared errors compared to the measured values for f(t). EQ. 3 1 is linear in Xp yp Zp ap bp and Cj 
so minimum norm methods, such as least squared errors, may be employed to solve for the 
unknowns. The number of frequency intervals and the lower and upper frequencies for each 
function segment are values that are selected, and thus, are known. Similarly the times when the 
measurements are taken are known. These values are inserted into the equation in FIG. 14 and 
the result is a set of coefficients that multiply each of the Xpyp Zp Op bp and Cp This set of 
coefficients forms a rectangular matrix ||A||. Designate the vector of measured values of f(t) as If]. 
The resulting matrix equation is EQ. 32. 

I|A|| |d| = |fl (32) 

Matrix solution methods known to those skilled in the art can be used to solve for |d|, 
which contains the values for Xp yp Zp Op bp and Cp One such method is using the Moore-Penrose 
generalized inverse of ||A||. Designating ||P|| as the Moore-Penrose generalized inverse of ||A|| 
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leads to EQ. 33. Singular value decomposition is another method that can be used. 

l|P|||fl = |d| (33) 

For the situation where the measurement times and frequency intervals are predetermined 
and reused multiple times only one matrix inverse needs to be calculated, which offers 
considerable economy of calculation. The matrix ||A|| is based on the predetermined [Lp Hj] and 
predetermined measurement times and is calculated in advance. The Moore-Penrose generalized 
inverse, ||P||, is also calculated in advance arid stored for later use, such as in the nonvolatile 
memory of a computation means such as using a hard disk or other magnetic media or on a CD 
ROM or other optical storage media, or nonvolatile RAM or other electronic storage media. 
After ||A|| and ||P|| have been calculated and stored the values for ||P|| can be retrieved and used to 
calculate |d| whenever such calculations are desired for a vector of measured values of f(t), [f] 
using EQ. 33, which leads immediately to values for the coefficients Xp Zp Op bp and Cp 

As described below, the formulation becomes substantially simpler when the lower and 
upper bounds for each interval, [Lp Hj] all have the same width, but they do not need to have the 
same width or be equally spaced. No restrictions are placed on the measurement times being 
equally spaced. The flexibility available for selecting the bounds for each interval and the 
measurement times allows selecting nonxmiform intervals, such as using Chebyshev points or 
expanded Chebyshev points. 

Methods other than the least sum of squared errors can be used to calculate the values of 
the unknown coefficients . For example, the least sum of absolute value of errors (LAV) or 
nonnegative least squares (NNLS) can be used. As is known to those skilled in the art, LAV can 
be implemented such that constraints are placed on the values of the variables, as is also the case 
with NNLS. Such constraints are particularly useful when linear function segments are used 
because these readily lend themselves to constraints that prevent one or more of SDC functions, 
such as G(f), from having any negative values. 

For the case where spectra are continuous, such as the Doppler spectra from flowing 
blood, the jth function segment needs to connect to function segment j-L Requiring continuity 
means that the endpoint of the (j-l)th function segment must be the starting point for the jth 
function segment, which leads to the lower bound of the jth segment being equal to the upper 
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bound of the (j'l)^h segment. Therefore, EQ. 34 holds for piecewise continuous functions. 
^y = ^-i (34) 

As demonstrated below, this restriction means that the value for c. is determined by the 

values of , , b._^^ c._ ^ and fi, . Therefore, c. is not an independent variable when the 
5 quadratic function segments are piecewise continuous. 

Spectra from processes such as flowing blood are smooth. Smooth transitions between 
quadratic function segments requires that the slope at the end of the (j'O^h segment match the 
slope of the jth segment. As will be shown below, this constraint removes the set of bj from 
being independent variables. 

10 The following discussion uses as examples the ^^y? , and the same discussion also 

applies to the ; Include an offset for the frequency at the beginning of the interval over 

which a function is being used. The equation for the jth function segment that goes through a 
specific point is derived as follows. Designate the specific point for the lowest frequency and the 

amplitude at the lowest frequency by using the ordered pair ( L., P . ). p. is the beginning 

15 amplitude for the jth function segment. The jth function segment has its own coefficients bp 
and intercept Cj , 

The starting point for the j/// function segment has the same coordinates as the endpoint 
for the {j'l)th function segment. Designate the endpoint for the {j-l)th function segment as 

( i» 1 ). ^y_ 1 is the ending amplitude for the (j-ljth function segment. The starting point 

20 for they//? function segment is the endpoint of the j-lth function segment. The following 
equation gives the amplitudes for the jth function segment of the SDC function when using 
piecewise continuous quadratic functions to represent the SDC function. 

G. = a.{J-H._^)\b.{f-H._^) + c. (35) 
However, c. =e ,so 

7 7-1. 
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G. = «.(/-//._,)' + 6.(/-//._,) + e._, (36) 

For the special case when j-1 the starting point has an amplitude of 0 at the frequency 

fif^ , where is defined earlier. Clearly, a recurrence exists that gives the . The basic 
equation is: 

e-a.iH-H. j\bXH.-H, J + e, ^ (37) 

This equation gives the logical result that the height of the endpoint of the jth function 
segment is the sum of the heights of all of the preceding endpoints plus the change in height that 
occurs during frequency interval j. Using the preceding equations to eliminate Cj from EQ. 35 
leads to EQ. 38. 

(38) 

G =a,{f^H. J +bXf-H. J+ y ((// -H J a + (// J ) 

J J J - ^ J J - ^ p p- \ y p p p - ^ p 

These equations that tie the endpoints together establish continuity between the piecewise 

functions. The piecewise functions are now piecewise continuous. The values for the H. are 
known because they are selected a priori, hi this set of equations the only free variables are the 

firstA/-7 b andthei\^ a . 

J j 

Constrain the slopes of the j-lth and the jth segments to match where the segments joint. 
This constraint leads to the starting slope of the jth function segment being equal to the parameter 

bj The b. for each segment depends only upon the values for the coefficients for the previous 
function segments. The end result is that each b. depends only on the previous values of the a.. 
The equation for b, is 



^=2 

J 
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(39) 
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For the first function segment the initial slope is zero because it starts tangent to the f- 
axis. Therefore, ^, = 0 . This situation becomes clear by examining the case where j=l . 

G^^a^U-nJ (40) 

Using EQ. 39 to eliminate b from EQ. 38 leads to EQ. 41 . This equation is the general 
expression for the jth function segment for piecewise continuous quadratic function segments. 



J 
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The endpoint of the final function segment is at frequency . The SD function needs to 
be zero at this point and the slope also needs to be zero. These requirements also require the 
SDC functions, such as G^, to have zero amplitude and slope at . Two boundary conditions 

exist at , but only one free variable, ^ , remains, therefore need to remove another degree 
of fi"eedom. Set the equations for the boundary conditions that equate the slope and amplitude at 
the end of the Nth segment equal to each other (they both equal 0) and it becomes clear 

that J is not a free variable, either. Therefore, when have N function segments only N-2 

values of |a| need to be found. 

The equations become considerably simpler is equally spaced function segments are used. 
Assume that each of the function segments covers a frequency interval of equal width. Select a 

value for the lowest basis firequency, H , and then index by up fi-om it for other fi*equencies 

in the basis. 

With equally spaced frequency intervals have 

^g-r^g-2=^ (41) 

And the general expression for the jth quadratic function segment becomes EQ. 42. 
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The unknown variables are the o^. Represent these variables as the vector |a|. When 
impose continuity and matching slopes at the knotpoints where the function segments join the 
values for the (N-l)th and Nth values of |a| stop being independent variables. The expressions for 
these variables are 

N-l 

«;v-,= Z«,(P-^) (43) 



(N-l 
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The equations that give ^ and confirm that the last two coefficients are 

completely specified by the first N'2 coefficients and, thus, are not independent variables. These 
values lead to the last two function segments working together to make the final amplitude of 
G(/) become zero and with zero slope with smooth transitions at the knotpoints of the final two 
function segments. 

Expressions for the (N-l)th and Nth function segments in terms of the first N-I elements 
of |a| are 

(45) . 
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X (n^ i-(2j- l)H^+2f-2H^)a.+ (-H^-(N-2) +/)' a. (j - N)) 
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Similar relationships exist for F(f), Fj and |d|, including 

(47) 



X(-(2/7-l)//5t2/-2//o)x, +(-H^-(j-l)H^+f) X. 



(48) 



l' (H,(-(2j-l)H^ + If- 2 H„) X + (-H^-(N-2)H^ x U-N)) 

N-2 

F^=j;,(-x.(-^+NH^+H,f(j-N+l)) (49) 

J = ^ 

Substituting the preceding equations into EQ. 22 and integrating and substituting to 
remove the N-lth and Nth Oj and Xj leads to EQ. 50, This equation is the final equation for the 
function that approximates g(t) when using equally spaced quadratic function segments that are 
piecewise continuous and constrained so that the slopes at the knotpoints are equal and that 
slopes and values at the beginning and end of the SD function are zero. EQ. 50 is a general 

(50) 



8(0 -E<0 



i^-2 
2 
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1 

g (2 ( W- J) « f (fl^ + (N- 2) /Tj)) + 2 cos(3 II < (/Jfl +J H^) + 2 ( 1 -y) cos(2 « K/^ 



formula for calculating g(t) that is linear in the variables x. and a. . 

The times when measurements are made are known, as are the starting basis frequency, 
Hq and the width of each frequency interval, H^. After these substitutions are made the result is a 

set of equations consisting of coefficients on the unknown N-2 unknown variables a. and the N-2 

unknown variables x. , The coefficients can be extracted and put in a rectangular matrix ||A||. ||A|| 

will have M rows and 2(N-2) columns where: 

M = number of measurements, with the measurements being i = 1 . . M 
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N = number of function segments (which is the same as the number of basis frequencies). 

Set ||A|| so that the first N-2 columns contain the coefficients for Gj and the last N-2 
columns contain the coefficients for the Xj, The equation for the first N-2 columns of ||A|| is 
EQ. 51. The equation for the last Nr2 columns of ||A|| is EQ. 52, These equations assume that 
the envelope function is known, which is often the case, such as if the interrogating function is a 

(51) 

^i,r\ (2 ^^^0 +-/'^6) ^/) - 2 cos(2 71 {H^ + 0 - 1 ) H^) t.) 
+ 2iN-j-l)cos{2nL(H^ + NH^)) 
+ 2{N-J)cos(2nt.(H^ + (N-2)H^)) 
+ (-4N+4j+2)cos(2nt.iH^+iN- l)H^)))/iK\L ) 

(52) 

A,j^^_, = lE(t.)(-2(N-j-\)smi2nt.(H^ + NH,)) 

-2sm(2K(H^+JH^)t.)+sm(2nt.(H^+iN-2)H^)){-2N+2j) 
+ (4 N- 4j-2) sin(2 n t.iH^+iN- I) H )) 
+ 2 sm(2 n iH^ + (j - I) H,)t.))/ (n't.) 

continuous sine or has other known characteristics. If the signal being analyzed has not been 
modulated by an interrogating function then set the envelope function E(t)=l and the method will 
decompose the complete signal and not attempt to remove the effects of the interrogating signal. 

Designate the values of the (currently unknown) N-2 values of as |a| and designate the 

(currently unknown) N'2 values ofx^ as |x|. Designate the vector |d| as being |a| stacked over |x|, 

thus |d| contains the unknown variables Ojj' - I , . N-2 followed by the unknown variables Xj, 
y = . . N'2. The values Idl can solved as described earlier, using, for example, the Moore-Penrose 
generalized inverse. 

FIG. 15 through FIG. 23 illustrate the method just described using an example. The 
example uses a small munber of measurements, 20, and seven quadratic fiinction segments. FIG. 
15 lists the parameters used to calculate the ||A|| matrix shown in FIG. 16. The Moore-Penrose 
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generalized inverse of ||A|| is shown in FIG. 17. Only ||P|| is needed to calculate the values of |d| 
from measured values of f(t). The 20 values of g(t) used in the example were calculated using 
the G(f) and F(f) shown in FIG. 19. A successful calculation will produce estimated G(f) and 
F(f) close to those shown in FIG. 19. FIG. 18 illustrates the envelope function used to generate 
the g(t). g(t) is illustrated in FIG. 20 along with the points where measurement values were 
taken. FIG. 21 illustrates what g(t) looks like after removing the effects of the envelope function. 
Contrary to the appearance of FIG. 21, this function has a period much longer than what is shown 
in the figure. The apparent periodicity of the signal at 2 MHz is caused by the narrow range of 
the spectrum that underlies g(t) and f(t). 

The values of |a| and |x| were calculated by multiplying ||P|| times the vector of measured 
values of g(t). The resulting estimated SDC functions, G(f) and F(f), are shown in FIG. 22 and 
FIG. 23 along with the actual G(f) arid F(f). Each figure shows only one line because the 
estimated values are so close to the actual function values. With the estimated SDC functions 
matching the actual SDC functions it is tautological that the estimated SD function will match 
the actual SD function. 

This example illustrates some key features of the subject invention that distinguish it 
fi-om other methodology inventions. (1) When using precalculated ||P|| the calculations of the 
spectra involve a small number of multiplications and additions, which means that the 
calculations are not memory intensive and can be fast. The calculations do use numbers that can 
be close to the same size to the precision of the calculations can require on the order of 1 5 or 
more digits of precision. (2) The example discriminated the spectrum using a relatively small 
number of measurements. (3) The size of ||P|| increases linearly with the number of 
measurements, therefore if use precalculated ||P|| the computational cost of using more 
measurements increases only linearly. (4) The size of ||P|| increases linearly with number of 
measurements, therefore if use precalculated ||P|| the computational cost of using more basis 
frequencies (function segments) increases linearly. (5) Amplitude, SD, and SDC functions are 
calculated directly using rational functions, including using rational real-valued functions and 
rational real-valued continuous functions. 
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The solution method illustrated using linear algebra has practical limits for the number of 
basis frequencies that can be used. The basis frequencies being used are, by design, not 
producing orthogonal basis functions. Therefore, ||A|| becomes increasingly singular as more 
basis frequencies are added, which tends to make the problem increasingly ill-posed when 
digitization errors exist in jf]. However, other variations in accordance with the present invention 
are described below that remove the problem of ||A|| becoming increasingly singular. 

When using xmiformly spaced piecewise continuous quadratic function segments the SDC 
functions can be calculated using EQ. 53. EQ. 53 uses Qj and will calculate G(f). Substituting Xj 
into EQ. 53 will lead to calculating F(f). EQ. 53 was used to calculate the distributions in FIG. 
22 and FIG. 23. 

(53) 



X Heaviside (/- //, - (y - 1 ) H^) Heaviside (//, +y if^ -/) 
X (-(2/7- l)//,+ 2/-2/f^)a^ 



Heaviside (f-H^-(N-l)H^) Heaviside {H^-f+NH^) 



The equations needed for using equally spaced piecewise continuous quadratic function 
segments are summarized in FIG. 24. As stated earlier, functions other than quadratic functions 
may be used. FIG. 25 summarizes the equations needed for using equally spaced piecewise 
continuous function segments. The functions listed in these figures are examples of using the 

15 general method of the present invention and demonstrate its general use. As stated earlier, the 

functions do not need to be equally spaced, continuous, or polynomials. The basis functions used 
in FIG. 4 are B-splines. Among the other functions that can be used are those containing 
trigonometric and exponential functions. 

Measured values are necessarily digitized to some precision. Digitization introduces 

20 errors because the actual values of the function are not known. Instead, the values available have 
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been moved to the centers of regions bounded by the measurement precision hmits. These 
regions are discrete regions into which measurements are placed during the digitization process. 
These regions into which measurements are placed are the digitization bins. This shifting of 
actual values to digitization bins is illustrated in FIG. 26 for the case where the number of 
digitization bins is four. The solid line is part of the plot of a function consisting of three discrete 
frequencies. The crosses are the actual values at selected measurement times. The circles are the 
digitized values that occur when have four digitization bins. SD functions are necessarily 
calculated from data that have been digitized. The digitized data are an approximation of the 
actual function values, yet the desired SD functions are those that produced the actual functions. 
The present invention includes a means for using digitized data to estimate actual SD functions. 

The method described above will generate SD functions that accurately estimate the 
measured values of f(t) that have been digitized. The reproduced values are the digitized values 
of f(t) and the SD functions are those bounded by the lower and upper basis frequencies. When 
using digitized data the estimated SD functions for the digitized f(t) become increasingly close to 
the spectra for the actual (nondigitized) f(t) as the precision of the digitization process increases. 
As the precision of the digitization process increases the measurement error decreases and the 
problem becomes less ill-posed. The method described above will accurately estimate SD 
functions for the special case where the frequencies of the function are known, even when the 
amount of digitization error is extreme. This ability is illustrated in FIG. 27 through FIG. 29. 
These three figures show, the actual function values and the estimated function values when 
coefficients are calculated by minimizing the sum of least squared errors, FIG. 27, and by 
minimizing the sum of least absolute value errors, FIG. 28. FIG. 29 shows the results of both 
least squared errors and least absolute value errors along with the actual and digitized function 
data. The estimated results shown were calculated using the digitized data.from the four 
digitization bins illustrated. 

Typically, the frequencies that compose the function of interest are not known and an 
approach different than minimizing some objective measure of goodness of fit between estimated 
and measured (digitized) values is needed. The digitization process causes the measured values 
of f(t) to be moved from the real (and unknown values) to values that correspond to the centers of 
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the bins bounded by the lower and upper thresholds specified by the number of bins into which 
data may be placed and the lowest and highest values of the total measurement range. The main 
consequence of the digitization process is that the measured values have almost always been 
shifted to values that are either higher of lower than the actual values. These shifts tend to 
5 introduce spurious frequencies into the measured data. The spectral decomposition methods that 
use objective measures to minimize differences between digitized values and estimated values 
calculates SD functions for the data presented to it (the digitized data), not the actual function of 
interest. Unless steps are taken in addition to those already presented, when the data have been 
sufficiently corrupted by the digitization process the resulting spectra no longer portray the 

10 spectra of the nondigitized data. 

An example of calculated G(f) and F(f) with digitization errors when no steps are taken to 
account for digitization are shown in FIG. 30 and FIG. 3 1 . Except that they are based on 
digitized rather than actual values of f(t), these figures are based on the same data and used the 
same ||P|| as FIG. 22 and FIG. 23. The spectra for the digitized data can be used to calculate 

15 estimated values for the measured digitized data and the estimated values are the same as the 
actual digitized values, confirming that the calculations divined the correct spectra for the 
digitized data, while also confirming that spectra of the digitized data do not represent the spectra 
of the underlying fiinction that was digitized. 

Dithering the digitized data will improve the results, but dithering alone is often 

20 inadequate. Other steps are often necessary when processing data that have been digitized. 

Increasing the number of digitization bins will improve the results. However cost and 
availability of suitable electronic components can limit the amount of such increases that are 
practical. Other steps besides increasing the number of digitization bins are often necessary. 
One aspect of the present invention involves estimating SD and SDC fimctions in the 

25 presence of digitization errors by using error absorbing functions. Digitization causes more 

abrupt changes in the measured values of f(t) than actually occur. These abrupt changes manifest 
themselves in the frequency domain as the introduction of fi"equencies that do not exist in the 
actual f(t). Error absorbing functions are fimctions that are active outside of the frequencies that 
exist for the actual f(t) and that account to at least some extent for the spurious fi-equencies added 
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by digitization. One example of an error absorbing function is including one or more SDC or SD 
functions to the basis and selecting the frequencies for the additional functions so that they are 
outside of the interval for the actual f(t). For example, the starting frequency for the absorbing 
functions can be a multiple, such as two more, of the highest basis frequency expected for the 
actual f(t). When the error absorbing function is set up using the same functional form as those 
used for the SDC or SD functions then the same objective measures may be used to calculate 
solutions, such as minimizing the sum of least squared errors or LAV. Error absorbing functions 
can be used in conjunction with other methods such as dithering. 

Another aspect of the present invention involves estimating SD and SDC functions in the 
presence of digitization errors by using SD or SDC functions that have one or more properties 
that can be adjusted using an objective measure, such as by minimizing or maximizing the 
properties. The properties being adjusted may have a useful geometric interpretation, but such 
interpretations are not required. The adjustment of the properties may be subject to one or more 
constraints. Following are some examples of this aspect of the present invention. This aspect is 
particularly useful when the number of digitization bins is small enough to produce erroneous 
results and dithering does not produce correct results. 

The SDC functions calculated using data significantly corrupted by digitization, as 
illustrated by FIG. 30 and FIG. 3 1, are known to be incorrect because G(f) is not always 
nonnegative. G(f) must be strictly nonnegative because it is based on the cosine of the phase 
angle, and the cosine function is an even function. Therefore, one aspect of the present invention 
is to calculate SD functions by constraining the calculated results so that one or more SDC 
functions or the SD function are strictly nonnegative. 

For example, when G(f) is modeled using uniformly spaced piecewise continuous 
quadratic segments EQ. 54 is the general expression for the constraints for the Qj that prevent 
G(f) from being negative anywhere. Only y=2 . ,.N-1 for a total of constraints are needed 
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because if these constraints are met then the first and last function segments are also forced to be 
nonnegative. EQ. 54 accounts for the minimum values of a function needing to be located within 
the interval for which the function segment is defined for its contribution to G(f). If G(f) is 
modeled as linear segments then the nonnegativity constraints for G(f) are expressed by EQ. 55. 

SD and SDC fiinctions have geometric properties. Geometric properties are properties 
that are affected by measures such as size, dimension, length, area, shape, and curvature. For 
example, the SDC functions depicted in FIG. 30 and FIG. 31 for the digitized f(t) have much 
larger amplitudes for the ones for the actual f(t). Similarly, the areas bounded between G(f) and 
F(f) and the frequency axis are much larger for the case with digitized data compared to 
nondigitized data. The bounded area of a spectral component function, G(f) or F(f), is the total 
area between the frequency axis and the curve above the frequency axis added to the absolute 
value of the area between the frequency axis and the curve below the axis. Bounded area is a 
geometric property because it is affected by the area and shape of the SDC function. The heights 
of peaks and valleys of functions and the number of such features are also geometric properties. 

Curvature is a geometric property. Curvature is any combination of one or more 
functions that are first and higher derivatives of SD, SDC, or amplitude functions. One type of 
curvature when using piecewise continuous quadratic function segments is the squared 
magnitudes of the vectors of the coefficients |a| and |x|. This measure of curvature is called the 
quadratic curvature and it is the sum of the squared second derivatives of the quadratic function 
^segments. 

Similarity is another geometric property. Similarity is any combination of one or more 
functions of the parameters that characterize two or more SDC functions such that the parameters 
of a first and second SDC function are related. For example, the squared differences between , 
parameters that define G(f) and F(f) is a similarity. If G(f) and F(f) are characterized using 
piecewise continuous quadratic function segments described earlier then the sum of squared 
differences between corresponding elements of |a| and |x| is another function using similarity. The 
geometric interpretation is that the similarity property relates to similar shapes, slopes, or 
changes in slopes between two or more SDC functions. 
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Another aspect of the present invention involves calculating one or more SDC functions 
or SD functions while adjusting one or more measures of geometric properties. For example, in 
accordance with this aspect of the present invention, the arc length of one or more of the SDC 
functions or the SD functions can be calculated so that the arc length measure meets one or more 
5 objective measures, such as being minimized. Similarly, in accordance with this aspect of the 
present invention, one or more SDC functions or SD functions can be calculated while adjusting 
one or more measures of the boimded area of one or more of the SDC functions or the SD 
function so that the area measure meets one or more objective measures, such as being 
minimized. Similarly, in accordance with this aspect of the present invention, one or more SDC 

10 functions or SDC functions can be calculated while adjusting one or more measures of the 

bounded area of one or more of the SD fimctions or the SD fxmction so that curvature meets one 
or more objective measures, such as being minimized. 

Calculating SDC or SD functions while adjusting geometric properties can be 
accomplished using various techniques. One technique is to minimize the measure of the 

15 geometric property , such as to minimize a measure of arc length, bounded area, or curvature of 
the function. When G(f) is based on evenly spaced continuous quadratic function segments one 
expression for a measure of the total are length is EQ. 56. The same equation calculates a 
measure of the total arc length for F(f) when Xj are substituted for a,. 
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20 EQ. 56 is based on the square of the total arc length and it has had constant terms, such as 

the frequency increment, removed from the right hand side. Therefore, it represents a 
measure of the total arc length of an SDC function, rather than being the total arc length. These 
forms are independent of the frequencies being used in the basis. For example, for the case 
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where N=5 the following expression is the measure of the total arc length of G(f). 

(56) 

4 2 : = 5 + ^2 + ^3) (-4 a^^3a^-2a^) 

+ 2 (-4 - 3 ^2 - 2 a^f + 4 (a, + ^2 + 0^)^ + 0^^-^ 3 (a^ + ^2) tf, + a^^ 

2 

+ 3 (aj + ^2 + ^3 ) + ^2) + ^'3 

EQ. 56 is minimized when its gradient is calculated with respect to each of the 
coefficients, the and set equal to zero. EQ. 56 is quadratic in the coefficients, therefore the 
gradient with respect to each coefficient will be linear in the aj. For example, when N=5 the 
gradient is EQ. 57. Identical expressions exist for F(f), except that they use the variable x instead 
of a. 

(57) 

[46 + 30 ^2+ 13 a^ - 0, 30 a, + 22 ^2 + 10 ^3 = 0, 13 10 + 6 a^ = 0] 

The gradient equations are linear, therefore, they can be included as additional rows in the 
linear regression matrices used earlier, or in other calculation methods, such as LAV. The values 
used in the additional rows of ||A|| are the Hessian of the measure of the arc length measure. For 
example, if N=5 then the additional rows for the columns that multiply the vector |a| are 
"46 30 13" 
30 22 10 (58) 
13 10 6 . 



The additional rows in the vector \f\ will have some (nonnegative) values based on the 
underlying structure of the G(f) or F(f), or the additional rows in |fl can be set equal to zero to 
calculate values for |a| and |x| that minimize the squared errors between the estimated \f\ and the 
measured |f| while simultaneously reducing the total arc length of G(f) and F(f). The number of 
gradient rows is limited by the number of quadratic function segments selected to model G(f) and 
F(f). The number of gradient rows will necessarily be much smaller than the number of rows 
based on measured data. The contribution from the gradient rows can be scaled by multiplying 
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each row by the a constant. 

The example above used arc length as the geometric property. Other geometric properties 
can be used, such as bounded area, curvature and similarity. The measures of these properties 
also produce quadratic expressions that yield Hessian matrices with constant terms. 

When SDC functions are calculated using least squared errors and rows for arc length 
measures or other geometric properties with quadratic measures are included and set equal to 
zero the resulting SDC functions possess the general shape of the nondigitized SDC functions, 
however the curves are riiuch flatter than the actual SDC functions. As more function segments 
are included the general outline of the actual functions tend to become more apparent, but the 
results remain approximate. These approximate results are quickly calculated because the ||P|| 
matrix can be calculated in advance and stored for later use. These approximate results are 
useful starting points for more refined calculations of SDC and SD functions. 

The values for parameters that characterize SDC and SD functions can be calculated from 
the gradient equations, such as the equations for the gradients of the arc length measures 
presented above. When the parameters, such as |a| and |x| when using equally spaced continuous 
quadratic function segments, are calculated using gradient equations the SDC functions tend to 
change systematically, which is beneficial when improved SDC and SD functions are being 
calculated from starting values of parameters. 

In general, the values for the gradient terms for G(f) and F(f) that would be used in |f| are 
not known and the result is estimated spectral component functions that have amplitudes that are 
too small and too flattened. However, the results can be used as starting points for additional 
steps that can be taken to improve the results. 

As described' earlier, other geometric properties of SDC and SD functions can be used 
besides measures of arc length. The bounded area of SDC and SD functions is one such 
geometric property. FIG. 91 is an expression for a measure of bounded area of G(f) when 
piecewise continuous equally spaced quadratic functions are used as the basis functions. These 
expressions are based on the squared area between G(f) and the frequency axis and are quadratic. 
EQ. 59 is an example when N=5. 
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(59) 



MAs , 3 = 251 + 158 + 33 + 673 + 738 a, + 217 



The gradient term is linear in terms of |a| and can be included in a least squared errors 
calculation as was done for the arc length measure, or in other calculation methods, such as LAV. 
The additional rows for the ||A|| matrix are the Hessian of the arc length measure. For example, 
when N=5 the Hessian is 

"1346 738 251" 
738 434 158 (60) 
251 158 66j 

Similar equations for F(f) exist for the measure of area expressions and are the same as 
the ones for G(f) except that a is replaced by x. 

Nonnegativity and measures of arc length and bounded area are two examples of 
10 characteristics of SDC and SD functions that can be used to constrain the results of calculations. 
Other characteristics besides these can be used, such as those that are based on slopes, number of 
modes, relationships between mode peaks or valleys, ciirvature, similarity and other geometric 
properties. 

The digitization process provides digitized values for f(t), but the actual values for f(t) 
15 can be anywhere within the digitization bin within which a digitized value falls. The probability 
distribution for the values within a bin is uniform: all values are equally likely to be correct. 
Consequently, no reason exists to prefer one value over any other. One aspect of the present 
invention is to calculate SDC and SD functions by constraining the calculated results to produce 
calculated estimates of f(t) that are within predetermined limits of digitized values. A 
20 particularly useful set of limits are the lower and upper limits for each of the digitization bins. 
This approach contrasts with standard least squares estimates, such as those of parametric 
methods, which minimize squared errors between actual and estimates, but do not constrain 
predicted results to be within predetermined boundaries based on the precision of the digitization 
process. Fourier analysis methods fit data to functions in a least squared errors sense, such as 
25 Fourier series fitting data to trigonometric polynomials in a least squares sense, in contrast to the 



45 



aspect of the present invention that requires results to be within limits, but does not bias the 
results to come close to any particular value within the limits. 

The measured values for f(t) have been digitized and the actual (nondigitized) values are 
not known. However, each of the digitized values for f(t) have known lower and upper limits. 
Each of the digitized values is a member of a digitization bin that has well defined lower and 
upper limits and any estimated value for a particular value of f(t) must be between the lower and 
upper constraint values. When all of the constraints for each of measured values of f(t) are 
considered collectively a set of digitization bin constraints exists. 

The digitization bin constraints and the nonnegativity constraints can be combined with 
the general objective of minimizing one or measures of geometric properties, such as arc length, 
bounded area, curvature or similarity. Assuming that equally spaced piecewise continuous 
quadratic function segments are being used to model the spectral component fimctions then the 
equations in FIG. 33, FIG. 34, and FIG. 35 can be combined mto objective functions to be 
minimized by using individual expressions or adding expressions together. All of the 
expressions in FIG. 33, FIG. 34, and FIG. 35 are convex in aj and Xp therefore a unique global 
minimum exists for each of them and the local minimum is also the global minimum. These 
expressions may be combined into a single objective function to be minimized, such as the 
expression for the objective function in FIG. 36. This expression is also convex in Oj and Xj 
because the sum of convex functions is also a convex function. 

Another aspect of the present invention involves calculating spectral densities from 
digitized data by minimizing or maximizing one or more functions that are continuous over at 
least part of the frequency domain. A further aspect involves, for such minimization or 
maximization, using one or more convex functions (when minimizing) or one or more concave 
functions (when maximizing) as objective functions in calculating SDC or SD functions. Using 
convex functions as objective functions for calculating SDC or SD functions is particularly 
beneficial when digitization errors exist in the measured values of f(t). Convex functions, unlike 
honconvex functions, have unique global minimum values that are also the local minima for the 
functions. These properties make finding the minimum values of convex functions much easier 
than finding the minimum values of nonconvex functions. Instead of minimizing functions that 
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define properties of G(f) or F(f) it is possible to maximize functions. When maximizing 
functions it is beneficial to use concave functions because they have unique global and local 
maxima that are the same. 

Another aspect of the present invention involves calculating SDC and SD functions by 
minimizing or maximizing one or more objective functions subject to constraints. Such 
constraints include using digitization bin constraints or nonnegativity constraints or using both 
digitization bin and nonnegativity constraints . Using digitization bin constraints is particularly 
advantageous. For example, if the objective is to minimize the arc length measure of G(f) and no 
digitization bin constraints exist then the shortest length is a straight line extending between the 
lowest and highest basis frequencies along the frequency axis. When digitization bin constraints 
are included then one or more of G(f) and F(f) must be nonzero if one or more of the values of 
f(t) are nonzero. 

Two digitization bin constraints exist for each measurement: each estimated value must 
be larger than the corresponding minimum value for the bin and the each estimated value must be 
no larger than the corresponding maximum value allows for the bin. As the number of 
measurements increases the possible shapes and dimensions that G(f) and F(f) can have become 
increasingly restricted. . 

Convergence to the G(f) and F(f) that approximate the actual G(f) and F(f) for the 
nondigitized f(t) is often facilitated by identifying the actual lower and upper frequency limits, fL 
and fH. The search for these limits may be started by starting with limits known to be wider than 
the actual limits and searching for the |a| and |x| that rneet the digitization bin constraints and that 
minimize the objective function. Such a search is often facilitated by included nonnegativity 
constraints that prevent G(f) from having values less than zero. 

The number of peaks for G(f) and F(f) may be know a priori, or they may be determined 
during initial calculatibns when nonnegativity constraints are used. If the number of peaks in 
G(f) are known to equal the number of peaks in F(f) then the number of peaks for F(f) can be 
determined using initial calculations. The method of the present invention quickly converges to 
solutions that correctly identify the number of peaks in G(f). When the nxunber of peaks in F(f) 
is known then constraints that force F(f) to have a set number of peaks can be included. Such 
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constraints can also be used to force G(f) to have a specified number of peaks. Similarly, 
constraints can be imposed to force a specified number of valleys to exist in G(f), F(f), or both of 
these functions. 

- In general, numerical methods must be employed to find the constrained minimum or . 
maximum values of the objectiveTunctions. The nonnegativity constraints can be formulated as 
convex or monotonic functions. The constraints on the number of peaks or valleys of G(f) and 
F(f) can be formulated as convex or monotonic functions. The digitization bin constraints can be 
formulated as convex or monotonic functions. The objective function can be formulated as a 
quadratic function, as was done above. Therefore, the problem to be solved is a convex 
programming problem. Many methods, including quadratic programming, penalty function, 
barrier function, and those based on Lagrangian methods, exist and are well known. Penalty and 
barrier function methods are particularly well suited to solving the problems associated with 
estimating SD functions using the methods of the present invention. 

The functional form of the envelope function E(t) in EQ. 23 may be known, although the 
values of the parameters may not .be known. For example, a Gaussian pulse may used (EQ. 2) 
but the exact values of its parameters may not be known. The parameters of the envelope 
function can be quantified as the parameters for the SD or SDC functions are quantified. 

The following description uses a specific medical device application for illustrative 
purposes. As will be appreciated, the present invention may also be readily employed in 
applications beyond the specific medical device application described, such use with living 
organisms other than humans or being used in applications other than medical. As will also be 
appreciated, the present invention may take on configurations other than those described and 
illustrated. For example, the means used for calculations may combine or separate functional 
units in configurations other than those described. As another example, the sequence of 
calculations may. be done differently than those described. As another example, the independent 
variable may be a quantity other than frequency. 

The general approach of the preferred embodiment consists of two steps: 
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(1) set up the problem and 

(2) solve the problem 

The present invention contains novel and useful means for both of these steps. 



The preferred embodiment minimizes a constrained objective function defined as: 

COF= + + NCV + PCV (61) 

Where 
• OF^ \d\\\Q\\\d\ 
Note: this is an algebraically equivalent form of the objective function comprised of the 
measures of the geometric properties of G(f) and F(f) 
where 

Q = Vi of the Hessian of the objective function, 

A. = the scale factor that makes ^^^^^ ^ — ^ have a suitable scale. ^ is calculated 

using the latest estimate for |d|. To start by forcing the solution toward meeting all constraints the 
contribution of OF to the value of COF can be reduced by using a value of ^ that makes 

kllleiiki 

have a value that is smaller than the magnitudes of CVV, NCV and PCV until the 

values of CVV, NCV, and PCV have started to diminish because constraints are being met. 
More preferred is using a value of ^ as being two orders of magnitude larger than 

\d\'\\Q\\\d\. 

A preferred implementation sets up the OF using a plurality of functions that when 
combined using one or more mathematical operations, including, for example, the operations of 
addition, subtraction, multiplication, and division, produces a combined function that spans a 
predefined frequency range such that one or more combined functions approximate the SDC 
functions. Even more preferred is for the functions comprising a combined function to be 
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polynomials. Still more preferred is to use real valued piecewise continuous functions. Yet more 
preferred is to use real valued piecewise continuous functions that are equally spaced along the 
frequency axis. Finally, even more preferred is to use piecewise continuous real valued quadratic 
functions that are evenly spaced along the frequency axis that have first derivatives with 
matching values at the points of intersection. 

For measuring the spectra from blood flowing in an artery digitization errors will be 
significant, therefore a preferred implementation solves for the SDC functions by setting up and 
solving a constrained optimization problem. Preferably, the objective function will minimize or 
maximize one or more functions that are measures of geometric properties of one or more SDC 
functions, G(f) and F(f), or functions that are derived from G(f) and F(f), such as the SD function 
and the phase angle. G(f) and F(f) are based on the cosine and sine of the phase angle, 
respectively, and the amplitude of the SD function. The preferred objective function (OF) will be 
convex. Even more preferred is that OF is quadratic in |a| and |x|. 

A preferred implementation sets up the problem so that it includes one or more 
digitization bin constraints and may include one or more nonnegativity constraints and may 
include one or more peak count constraints for G(f) or F(f). A preferred implementation sets up 
the constraints as penalty functions or barrier functions. 

Particularly useful for the present invention is to formulate the digitization bin constraints 
as penalty functions that have small values when constraints are met and ever increasing values 
as deviations from the constraint limits increase. The digitization bin constraints require that the 
estimated values of f(t) fall within both lower and upper limits, therefore the digitization bin 
constraint value function has two sides that slope down toward the digitization bin boimdaries. 
Such functions are double sided constraint boundary functions. Of even greater use are doubled 
sided constraint boundary functions that have slopes that can be adjusted as the search for the 
solution progresses. One such double-sided constraint boundary function is 




(62) 
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where x=deviation of estimated value for f(x) from the midpoint of the digitization bin 
p is a positive, even, integer 

As p increases the digitization bin constraint violation value function becomes steeper, as 
illustrated in FIG. 37. For the preferred embodiment the value of the exponent p increases during 
the solution search process from a small value (e.g., p=2) to a large value (e.g., p=large value). 
As p increases the digitization bin constraint value function starts with sides witii small slopes 
(e.g., p=2) and the slopes increase as then p increases until b(x) has a very small value at the 
constraint boundary because the sides have steep slopes (e.g., p has a large value). Typically, the 
final value of p will be large: much greater than 100, and better yet have values larger than 
1,000. Still better, is to have final values of p approximately equal to or larger than 10,000. This 

function, in the limit of /? = ~ , is a step function at the constraint boundary, w. 

The digitization process causes the actual value of |fl to almost certainly differ from the 
true value, but the digitized value is at the center of if the corresponding digitization bin and this 
information leads to M additional constraints. Extending the preferred CVV function to each of 
the measurements is now described. Let 

M = the number of measurements. 

|f] is the vector of measured values. |f) is a column vector with M values. 
Let CVV = digitization bin constraint violation value. The preferred embodiment of the 
expression to calculate this value is: 



M 



cFF=2;B,M,rki 



(63) 



/= 1 



where ' 




This expression is a double sided constraint violation value function. 



I A. I is the ith row of the matrix ||A||. 



T T 

\A.\ \d\ is the estimated value of the ith measiirement and \A.\ \d\-f. is 
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the difference between the estimated value and the measured value of the zth measurement. 

p is an even, positive, integer. If p=2 then there is a least squares fit. As p 
increases the constraint violation value function becomes steeper and its corners move closer to 

the digitization bin constraint values of f^ - ^ and /, y • 

The preferred embodiment specifies two ways to calculate NCV using one of two similar, 
but still different, fiinctions. Refer to FIGS. 38 and 39. These functions are based on the general 
function form of 

l^CV^p^GjJ.JfgJgJg^ (64) 
where 

fgy^ = gate function that tends toward 1 when the extreme value for G. is on the 

frequency interval where G, is active and otherwise tends toward 0. 

fg^ = gate fiinction that tends toward 1 when G.{f.^^) < 0 and that otherwise tends 
toward 0. 

> 

fg^ = another gate fiinction that tends toward 1 when Gjif. ^f)< 0 and that otherwise 
tends toward 0. 

f = fj^ , the normalized frequency at which the extreme value for the yth function 

segment occurs. The extreme value for Gj occurs when f is between -1 and 0. 

The preferred embodiment for fg^ is based on setting up two fiinctions that are 
multiplied together. 

1 ) When fj^ <0: Set up a fiinction that is 1 when the argument is less than 0 and 
otherwise is 0. 

2) When f,^ >-l : Set up a fiinction that is I when the argument is greater than -1 and 
otherwise is 0. 
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2((l+/)^+F)^ ' 
The product of the two functions covers -^<fj g^, <0. 



p-2 

J 



f2 



(67) 



This is the expression that makes the noimegativity constraint violation active when f. 

is on the frequency interval where Gj{f) is active. Refer to FIG. 40 for an illustration of this 

gate function using two tapered Heaviside functions. Tapered Heaviside functions are described 
in more detail below. 

The terms provide three features, k is a positive value between 0 and 1, and typically 

will be close to 1 . The exponent p is a positive integer. Therefore, is always non-zero, f 

of 1+/^ could be 0 so adding prevents a singularity from occurring. Asp increases k^ 

will decrease toward 0, but always be a positive value. / is a normalized frequency that will 

have values on the order of 1 . The presence of k^ at low values of p makes the gate functions 
rise slowly and have slopes that become steeper as p increases . At large p the gate functions 
approximate Heaviside functions, but without the singularity that occurs with the Heaviside 
function. These functions that approximate Heaviside functions, but that are tapered, are called 
tapered Heaviside functions. The tapering of these functions makes them much more useful than 
using Heaviside functions when searching for solutions. Using tapered Heaviside functions in 
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constraints is a feature of the present invention. The ability to adjust the slope of the gate 
function facilitates convergence during the search process. The process can start with low values 

of p and then increase p to force the solution to ever more tightly enforce constraints. The 

term provides a second useful feature besides providing a means of altering the slope of the gate 

functions, the presence of causes non-zero gradient and Hessians to occur, which are 

necessary for using Newton's method to find solutions. 

The preferred embodiment for fg^ is based on setting up one function. The second gate 
function is to have a value that tends to 1 when Gif. ) < 0 and that otherwise tends to 0. 
Following the same function form as just used, set up a function that is 1 when G if ) <0 and 

Otherwise is 0. Base the expression on / =j which will have a value of 1 when 

Gjifj ) <0 and -1 when G.if. ) >0. This base expression has a total amplitude of 2. The 
desired range for the function is [0, 1], so divide by 2 to obtain the correct span for the range and 

then add y to obtain a function that has the desired range. To avoid a singularity when 

) =0 and to have the expression to produce increasingly sharp results as change a 

parameter, include a constant in the denominator. Unlike the functions based on /. , , this 
function does not have a natural scale of 1. Therefore, include a scale factor, . The extreme 
amplitudes of each G.{f) are linear in the a. , and the other term in the denominator is 

2 

Gjifj ) , therefore a reasonable scale factor to use is the mean squared value of the aj in |a|. 

Another approach would be scale the by normalizing the function for which the spectral 
density function is being estimated. In that case, the natural scale factor would become 1 . For 
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the general case, the preferred embodiment uses the following expression for the second gate 



function: 



fg2=- 



1 G 



1 




(68) 



where 

0<k<l ^ > 

\<p 

N-2 

|Li : scale factor, \i == jq _2 
and 

G = ) , the extreme value for the yth function segment, whether or not it occurs 

over the region where Gjif) is active. 

Refer to FIG. 41 for an example illustrating the second gate function. This figure shows 
the tapered Heaviside gate function for the case when the extreme value of G(f) is less than zero. 

The preferred embodiment for fg^ is based on setting up one function. As was the case 
with the second gate function, the third gate function is to have a value that tends to 1 when 

Gj{fj ^^t)<^ and that otherwise tends to 0. Following the same function form as just used, set 
up a function that is 1 when G .(/ ) is in a valley and that otherwise is 0. The second 
derivative of the yth segment of G(f) is 2 a. . Therefore, when a. >0 the extreme is a valley. 
Using the same function form as just used: Set up a function that is 1 when a. >0 and 



otherwise is 0. Base the 



a. 
J 



which will have a value of 1 when a, >0 and -1 
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when a. <0. The base expression has a total amplitude of 2 and the desired range is [0, 1], so 

divide by 2 to obtain the correct amplitude and then add y to shift the zero. To avoid a 

singularity when a. =0 and to have the expression produce increasingly sharp results as a change 

. in parameter, include a constant in the denominator. Unlike the functions based on fj , this 
function does not have a natural scale of 1. Therefore, include a scale factor, ^ . The other term 

in the denominator is Qj , therefore a reasonable scale factor to use is the mean squared values 

of the a J in|a|. Another approach would be scale the Oj by normalizing the fiinctidn for which 
the spectral density function is being estimated. In that case, the natural scale factor would 
become 1 . For a preferred implementation use the following. This expression makes the 

nonnegativity constraint violation active only when Gjif.^^) is less than 0. It has a range of [0, 
1], the same as the expressions derived earlier. 



^^^"2 / 2 ' .„ "^2 (69) 



where 
0<k<l 



N-2 

2 



H : scale factor, |ii = ^ 



N-2 

The general expression for NCV, EQ. 64, has a leading weighting factor with the 

exponent x. When x=2 the product of the gate functions is weighted by Gj(fj , which 
produces a nonnegativity constraint value function referred to as NCVl. FIG. 38 presents NCVl 
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for the preferred embodiment case of equally spaced piecewise continuous quadratic functions. 
NCVl forces valleys of G(f) up to the frequency axis, and then stops influencing the position of 
the valley. If other components of COF are attempting to make a valley negative then an 
equilibrium position with the valley resting on the frequency axis results. This situation is 
illustrated in FIG. 42. When x=l in the general expression for NCV then the nonnegativity 
constraint value function referred to as NCV2. FIG. 39 presents NCV2 for a preferred 
implementation case of equally spaced piecewise continuous quadratic functions. 

One aspect of the present invention relates to the ability to determine the values for the 
lower and upper bound frequency limits, fL, fn- When these limits are not known, a preferred 
implementation uses a search to find them and this search uses NCVl as the nonnegativity 
constraint in conjunction with the quadratic curvatures of both G(f) and F(f) as the objective 
function being minimized and the digitization bin constraints. The search starts using a wide 
frequency interval. The results are then used to identify one or more narrower frequency ranges 
to be used as the intervals for another sequence of calculations to determine function parameters. 
This iterative process may be used as many times as needed until no further improvements occur. 
The search is often adequate after fewer than 10 iterations, and not using more than this many 
iterations is preferred because the successive improvements in the results decrease. This iterative 
search for the frequency limits is preferred for use when such limits are not known before 
calculations begin to find the best estimates for the SDC functions. After the frequency limits are 
determined then the NCV2 nonnegativity constraint is used when making the best estimates of 
the SDC functions. 

The peak count constraint value function (PCV), FIG. 43, is very similar to the 

nonnegativity functions. It uses two gate functions: one that tends to 1 when the a, <0 
(meaning that the extreme value is a peak) and another gate function that tends to 1 when the 

extreme value for G, is on the frequency interval where G. is active. Tapered Heaviside 

functions are used in the expression. This function essentially counts the number of peaks and 
squares the difference between the number of peaks and target number of peaks. The magnitude 
of the result is unity, the same as the magnitudes of the other constraints. As the PCV constraint 
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is met its value will decrease as the square of the difference. For these reasons, do not include 
any additional weighting factor. To force the term used in the NCV function more quickly to 

a small value, the preferred embodiment uses ^^'^''^ in the PCV function. 

A preferred process for calculating the solution that minimizes the COF is to use a 
multidimensional Newton's method search with backstepping. The basic sequence of the 
preferred sequence of calculations to estimate SD or SDC functions is listed in FIG. 44 through 
FIG. 46. This sequence includes the steps used to set up the calculation sequence and then to 
execute the calculation sequence. 

The ilustrated process demonstrates that a scale factor, A, is used in the calculations, but 
its value is easily calculated, unlike when scale factors are used with regularization, such as 
Tikhonov regularization. This process also illustrates that there is no need to estimate the order 
of a time series model, as is the case with parametric spectral analysis methods. The illustrated 
process further illustrates that the present invention does not require fitting data using a norm 
(such as least squared errors) that minimizes errors between estimated values and measured 
values, and that the method explicitly allows estimated values to take on only values that are 
within the precision of the estimated values. Explicitly allowing estimated values to take on any 
value within the corresponding digitization bins, but constraining solutions to only be selected 
that make estimated values be within the correct digitization bins, provides a better solution than 
other spectral analysis methods. 

The basic sequence of calculations in FIG. 44 through FIG. 46 requires starting with an 
assumed |d|. The starting values for the elements of |d| are often approximately known, in which 
case these approximate values should be used. In the absence of approximate values, the values 
in |d| should start with very small nonzero values to reduce the number of iterations needed for 
the solution to be found. In the absence of an approximate estimate for the values of |d|, the 
preferred starting values for |d| are less than 10"\ and more preferably less than 10"^, and even 
more preferably less than 10'^^, and yet more preferably less than 10"*^. 

A preferred objective function to use is the sums of the squared values of |a| and of |x|. 
This objective function is the same as the squared magnitude of |d|. This function will be 
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minimized. 

The basic sequence of calculations in FIG. 44 through FIG. 46 may be used with the 
complete set of measurements from the beginning of the calculations. Using a subset of the total 
measurements accelerates the calculations when only a poor estimate is available for the starting 
values of Idj. The method in accordance with the present invention converges quickly to 
reasonable values of |d| that define the SD or SDC functions, even when a small number of 
measurements is used. Starting with a very approximate estimate of |d| will usually quickly 
converge to a reasonable estimate, that can then be used as the starting |d| when the fiill set of 
rneasurements is used. When the number of measurements exceeds a modest number, for 
example exceeding 128, a preferred implementation starts with a subset of the measurements that 
spans the total measurement time interval. A preferred number of measurements for the initial 
calculations is fewer than about 4,096 with an even more preferred number being less than about 
512. A preferred number of measurements for initial calculations exceeds 10. A more preferred 
number of measurements for initial calculations is iii the range of 32 to 128. For example, 64 
measurements could be selected for initial calculations, 

A preferred implementation includes backstepping. This backstepping differs from the 
typical backstepping employed in Newton and quasi-Newton methods. The preferred 
backstepping reduces each successive trial change in |d| by an order of magnitude. The preferred 
approach is to try between 5 and 15 orders of magnitude reduction in steps of one order of 
magnitude in attempts to reduce the COF. Calculating the starting step is done in the usual 
manner with Newton's method, as is known. This calculation is computationally more expensive 
than attempting successive reductions in the step size as attempts are made to reduce the COF. A 
variety of refinements can be included after a step size has been found that reduces the COF. For 
example, the interval between the last two steps can be further subdivided to attempt to find an 
even lower COF. Such refinements can further include using a quadratic or other minimization 
to estimiate the step size that gives the lowest COF. 

In a preferred implementation the value of the exponent p present in tapered Heaviside 
functions and in the double sided constraint functions starts at small values and then is increased. 
All values of p must be positive even integers for it to be compatible with all of the constraint 
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value functions. A preferred starting value is p=2. In a preferred implementationt the final value 
of p depends on the desired sharpness of the double sided constraint function for CVV. A 
preferred implementation for calculating the final value of p follows. Define e to be the 
allowed deviation from 0 at the digitization constraint boundary. It is a value that effectively 
5 becomes the same as 0. When a constraint is less than or equal to e then the constraint is met. 
In a preferred implementation the value of e is much less than 1 and even more preferably the 
value of e is less than 0.001 and still more preferably the value of e is less than or equal to 
0.00001. Define a to be the fraction of the digitization bin width at which e is to occur, oc is 
how close to the edge of a digitization bin e is to occur. In a preferred implementation the value 
10 of oc is between 0 and 1. Even more preferred is the value of cx being between about 0.9 and 1. 
Still more preferred is the value of a being between about 0.99 and 1. Yet more preferred is the 
value of cc being between about 0.999 and 1. Calculate the value of the final value of the 
exponent p based on the fraction of the width of the digitization bin width to where deviation 
from 0 is e . The value of the final p is the even, positive, value that most closely satisfies the 



15 following equality ^ > ^^ich equals Therefore, a preferred 

implementation selects the final value of/? as the even, positive integer that most closely makes 



the equality p = , which is the same as EQ. 70. 



loglO(e) 



For example, if e = 0.00001 and a = 0.999 then the final value of p is 

20 pValStop := 11508 . 

The value of p is increased after each Newton method search with backstepping. A 
preferred increase in is a value that is the next even positive integer that is no larger than 10 
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times the value of p most recently used. More preferred is increasing jt; to a value that is the next 
even positive integer that is no larger than twice the value of p most recently used. Even more 
preferred is increasing p to the next even positive integer that is no larger than 3/2 the value oip 
most recently used. In a preferred implementation the minimum increase in p is calculated by 
adding 2 to the most recently used value of p. 

For each value of /? a Newton method search with backstepping is done until no 
significant changes in COF occur or until no significant changes in any of the elements of |d| are 
made or until the Newton search procedure has executed a specific number of search steps. In 
the preferred embodiment significant changes in COF or |d| that lead to further Newton search 
iterations are fractional changes that exceed 10'^ of the starting values and even more preferably 
fractional changes that exceed 10"^ and still more preferably changes that exceed 10'^. In the 
preferred embodiment the maximum number of Newton search iterations for a particular value of 
p is between 5 and 100 and even more preferably between 10 and 50 and even more preferably 
between 10 and 25. 

In a preferred implementation the number of basis frequencies, N, is between about 3 and 
about 50, and more preferably between about 6 and about 30. In a preferred implementation 
when the number of basis frequencies is greater than about 10 the calculations start with a value 
of N that is smaller than the final N and the SDC or SD functions for that |d| that are used to 
determine the starting |d| for the final N, Even more preferably, the starting value for would be 
a factor of the final N, For example, if the |d| is desired for N~21 and G(f) and F(f) are being 
approximated by quadratic polynomials, then the first |d| calculated would be for N=7, The G(f) 
and F(f) from the |d| from N=7 is then used to calculate |d| for N=21 by fitting 21 segment 
piecewise continuous quadratic segments to the |d| calculated using Such curve fitting can 
be done using standard linear regression techniques. This |d| is then used as the starting |d| for the 
set of calculations that use N=21 . 

A preferred implementation performs the calculations using an augmented computing 
system, wherein such augmented computing system employs one or more electronic, magnetic, 
or optical components to perform digital or nxxmerical calculations, and store or display the 
results of calculations. Augmented computing systems include devices that use one or more 
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computers, microprocessors, or microcontrollers configured with nonvolatile or volatile data 
storage or retrieval means such as magnetic, optical, or electronic data storage or retrieval means 
and further configured to display or save calculation results or pass the results to a further 
processing step or system. 

A preferred implementation saves the calculated results for use by other algorithms that 
rely on spectral density data, such as those that use such information for calculating flow rates or 
for diagnosing diseases. An example of such an algorithm and associated structure is disclosed 
in PCT Application No. WO 03/032556 , published April 17, 2003, which is incorporated herein 
by reference. Such storage can be done using methods that employ one or more of electronic, 
magnetic, optical, or hard copy implementations. 

Another version of a preferred implementation displays data using optical means such as 
numeric or graphical displays such as those using one or more electronic means or hard copy. 

FIG. 47 through FIG. 53 are example displays for calculations using the above procedure. 
All of these figures were calculated using 18 digitization bits, which corrupts the measured f(t) to 
the point where the spectral density estimates calculated using ||P|| are incorrect by more than an 
order of magnitude and G(f) has significant nonnegativity constraint violations. FIG. 47 through 
FIG. 49 were calculated using 50 measurements. The results are close to the actual SDC and SD 
functions using a small number of measurements, demonstrating that the method produces useful 
results with an economy of measurements. FIG. 50 through FIG. 53 were calculated using 200 
measurements. The results are closer to the actual SDC and SD functions, as they should be 
because the results are more constrained by the larger number of digitization bin constraints. 
FIG. 53 is the cumulative distribution calculated from the SD function. Cumulative distributions 
can be more useful than density functions for some purposes, such as calculating dimensions of 
tissue structures. One aspect of the present invention relates to its ability to generate continuous 
spectral density functions that are more readily converted into accurate representations of the 
cumulative spectral functions than methods that produce values at specific frequencies, such as 
Fourier methods. 

FIG. 54 illustrates an example of a medical application in accordance with the present 
invention. In particular, the illustrated system uses an ultrasonic transducer 1 to introduce an 
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interrogating signal 2 into a patient 3, such as for calculating flow rates or dimension related 
information, using a transform utility in accordance with the present invention. The transducer 
may be any of various conventional ultrasound units and may provide the signal 2 in the form of 
a pulse train. In accordance with the present invention, the transducer 1 may be disposed external 
to the patient 3, for example, in the suprastemal notch and may be directed so that the signal 
targets the ascending aorta. Depthwise targeting can be accomplished by processing return signal 
components received during the time window corresponding to the desired depth. 

The interrogating signal 2 backscatters from patient tissues 4 and a return signal 5 is 
detected by the receiving transducer 6. The receiving transducer 6 converts the mechanical 
energy (not shown) of the return signal 5 into oscillating electrical signals (not shown) that are 
introduced into the analog amplification and processing module 7. This module may perform a 
variety of functions including analog digital conversion, amplification, filtering, conditioning and 
other signal enhancement. The output from the analog amplification and processing module 7 is 
a digitized stream of data that is the input to the digital processing module 8. The digital 
processing module 8 converts the digitized signal data into frequency spectra, in accordance with 
the present invention as described above, and then into velocity data which becomes velocity 
vectors that the digital processing module 8 stores as vectors of digital dimensioned velocity 
profile data in a memory location (not shown) of the data processing module 9. The data 
processing module 9 converts the velocity profile data (not shown) into one or more data vectors 
(not shown) that represent a piecewise continuous dimensioned velocity spectral density 
function. 

The digital processing module 8 may contain one or more microprocessors and one or 
more data storage structures that enable it to make the calculations needed to determine 
Volumetric Flow Rates (VFRs) and derived parameters as decribed in the above-referenced PCT 
application. The data storage structures contain values for parameters used in look-up tables that 
facilitate calculations that avoid or reduce the need to evaluate special functions. Additionally, 
certain values such as a channel dimension, average flow rate or the like may be predetermined 
and stored in cache or other storage for combination with values obtained in connection with a 
later measurement process. FIG. 54 also illustrates that the system has a control module 10 that 
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receive$ setup parameters (not shown) from the user 11. Setup parameters can include data such 
as patient parameters, e.g., hematocrit, and on/off signals regarding when data reading and 
processing should begin and end. The control module 10 receives signals from the data 
processing module 9 regarding the status of the received data (not shown) and changes in the 
interrogating signal 2 that are needed to improve or maintain quality measurement data. The 
control module 10 optionally receives signals of biological activity 12, such as electrical signals 
related to cardiac activity. The biological signals 12 may be used to synchronize interrogating 
signals 2 to biological functions or to correlate calculated values such as VFRs to such biological 
functions. The control module 10 signals the transducer power module 13 when the interrogating 
signal 2 should begin and end. The transducer power module 13 provides the drive signal 14 to 
the ultrasonic transducer 1. The data processing module 9 generates output data 15 that go to the 
output module 16. 

The output module 16 generates visible displays 17 and audio devices 18 to inform the 
user of the system's operating status and the results of the system's interrogation of the patient 3. 

The output module 1 6 may also include ports for providing data to other instruments 
(such as an EKG) or processing systems, for example, via a LAN or WAN. It will be appreciated 
that the various processing modules described above may be embodied in one or more 
computers, located locally or interfaced via a network, configured with appropriate logic to 
execute the associated algorithms. Additionally, certain signal drive and processing components 
may be incorporated into the interrogating signal instrument such as an ultrasound system. 

The preceding description is not restrictive and presents an example of one 
implementation of some aspects of the present invention. 
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